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SECTION  I 


INTRODUCTION 

This  is  the  interim  report  for  work  on  the  project  Dynamic  Simulation 
of  Airborne  High  Power  Systems  for  the  period  of  August  15,  1979  through 
August  30,  1980. 

The  work  has  included  modeling  of  a  three  phase  ac  generator  including 
effects  of  speed  variation  and  field  saturation  in  sufficient  detail  to  accu¬ 
rately  simulate  start  up  and  severe  faults.  Data  are  included  for  a  specific 
machine  and  results  are  included  for  simulation  of  a  severe  fault  condition. 

Preliminary  work  on  modeling  a  three  phase  three-leg  core  type  trans¬ 
former  is  included  here.  Circuit  equations,  magnetic  core  characteristics 
and  the  electric  circuit  relation  to  magnetic  core  characteristics  under 
saturation  conditions  are  developed. 

The  subcontractors  work  on  the  modeling  of  the  ac  resonant  charging 
circuit  is  reported  in  Section  IV.  This  work  consists  mainly  of  determi¬ 
ning  an  SCR  model  which  adequately  represents  turn-on  and  turn-off  phenomena 
in  a  high  power  reasonant  circuit  application. 

The  period's  work  also  addresses  the  problems  of  system  variables  in 
nonlinear  simulation,  stiff  differential  equations,  and  system  simulation 
using  SPICE  2  and  SCEPTRE.  The  subcontractor  shares  responsibility  and 
work  in  this  area. 

The  sections  of  the  report  correspond  to  the  tasks  outlined  in  the 
project  work  statement.  Detailed  development  and  data  are  included  in  the 
appendices. 

The  results  of  a  generator  simulation  run  are  included  in  Section  VII. 

The  results  of  the  ac  charging  simulations  are  included  in  the  Section  IV. 

In  the  ac  charging  circuit  simulation  a  single  loop  charging  and  a 
double  loop  charging  were  simulated  using  SPICE  2.  A  turn-off  problem  was 
corrected  using  a  snubber  circuit,  a  premature  turn-on  problem  is  being 
investigated.  The  SCR  model  is  being  converted  to  SCEPTRE. 

The  subcontractor  is  Virginia  Polytechnic  Institute  and  State  University 
with  Dr.  Fred  C.  Lee  as  Principal  Investigator. 
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SECTION  II 


THE  GENERATOR  MODEL 


2.1  OBJECTIVE 

This  part  of  the  project  was  to  develop  a  mathematical  and  computer 
model  of  a  three  phase  ac  generator.  The  model  is  to  be  adequate  for  in¬ 
cluding  effects  of  unbalanced  transient  loads,  startup  transient,  speed 
variation,  field  excitation  variation  and  saturation  of  the  field. 

Most  of  the  models  existing  in  the  literature  (References  1,5,6)  were 
developed  for  representing  generators  in  power  system  stability  studies. 

In  power  system  applications  the  generator  is  usually  connected  to  a  large 
system  and  the  generator  is  constrained  to  operation  in  a  narrow  range  or  it 
is  removed  from  the  system  by  protective  devices.  In  the  airborne  high  power 
system  one  has  a  single  generator  supplying  an  isolated  load  and  the  machine 
may  be  subject  to  more  severe  transients  of  a  wider  range. 

Further,  in  a  power  system  the  concerns  of  the  modeling  for  stability 
studies  are  sustained  overload  currents  and  the  synchronising  of  the  genera¬ 
tor  with  the  system.  In  this  airborne  high  power  system  the  operation  is 
asynchronous  and  there  is  need  to  consider  short  duration  transient  pulses 
that  may  affect  the  operation  of  the  electronics  in  the  associated  load. 

In  order  to  have  a  model  sensitive  to  these  concerns,  more  detail  is 
included  than  is  usual  in  power  system  stability  studies.  The  model  used 
has  these  additional  features:  direct  and  quadrature  damping  effects; 
variable  speed;  saturation  effects  including  both  variation  in  inductance 
and  variation  in  3L/3i;  more  accurate  developed  torque  formulation;  formula¬ 
tion  in  a  form  such  that  the  prime  mover  model  and  the  field  control  model 
may  be  added  if  these  become  available  at  some  later  time. 

2.2  GENERATOR  MATHEMATICAL  MODEL 

Appendix  A  shows  the  generator  equivalent  circuit  and  establishes  the 
notation  and  time  and  space  references  used  in  the  modeling. 


In  vector  and  matrix  form  the  equation  for  the  circuit  model  is  given 
either  by  Equation  (43)  or  (46)  and  repeated  here  as  Equations  (1)  or  (2). 


»■«♦&  (1) 

V  =  RI  +  £  (LI)  (2) 


The  torque  equations  relating  electrical  to  mechanical 
given  by  Equations  (55)  and  (56)  and  here  as  Equations  (3), 


variables  is 
(4),  and  (5). 


d2e 

dt2 


=  j[WB  af-] 


(3) 


where  the  Tg  terms  represents  the  so  called  electrical  torque  and  corresponds 
to  the  power  converted  from  mechanical  form  to  electrical  form,  Pg 


Concordia  (Reference  7)  shows  that  the  Tg  is  defined  by 


Te  =  \  [I]T  [&]  [I] 


(4) 

(5) 


Equations  (1)  through  (5)  basically  define  the  generator  model.  These 
equations  are  standard  forms  in  the  literature  (References  [1 ,3, 5, 6, 7]).  The 
model  used  in  this  work  differs  from  others  in  the  choice  of  system  variables 
and  in  the  details  of  defining  the  inductance  coefficients.  These  points  are 
discussed  in  Appendix  A  and  in  subsequent  sections  of  this  report. 


2.3  CHOICE  OF  GENERATOR  VARIABLES 

Choices  of  generator  state  variables  were  made  in  two  areas.  First,  the 
direct  phase  variables  were  chosen  over  the  more  traditional  direct  and  quadra¬ 
ture  axis  variables.  Secondly,  the  current  variables  were  chosen  over  the 
flux  linkage  variables. 

The  direct  phase  variables  were  chosen  over  the  traditional  direct  and 
quadrature  variables  for  three  reasons.  First,  while  Parks  (Reference  1) 
transformation  does  significantly  simplify  the  generator  model  equations  and 
hence  aid  in  their  solution,  it  is  necessary  at  each  numerical  integration 
step  to  transform  the  variables  back  to  direct-phase  variables  for  compatability 
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with  the  external  load  made  up  of  electronic  components  and  resonant  charging 
elements.  Any  expected  gain  would  be  greatly  diminished  by  this  transforma¬ 
tion. 

Secondly,  additional  logic  would  need  developing  to  represent  the  effect 
of  direct-phase  open  and  short  circuits  on  d-q  variables. 

Thirdly,  the  application  of  Parks  transformation  and  the  resulting 
simplified  equations  depend  on  certain  simplifying  assumptions  which  don't 
really  apply  if  saturation  effects  and  3L/3i  effects  are  significant. 

The  choice  of  the  current  variables  over  the  flux  linkage  variables  was 
dictated  not  by  generator  model  considerations  but  by  limitations  imposed  by 
models  of  the  electronic  components  of  the  load  (resonant  charging  circuit). 

The  developments  in  Section  V  and  VI  of  this  report  indicate  that  for 
appropriate  numerical  integration  methods,  use  of  the  A  variable  gives  less 
propagated  error  for  non-linear  inductances.  Further,  Nakra  (Reference  8) 
and  Manly  (Reference  9)  cite  inaccuracies  that  may  arise  due  to  noise  in  the 
317 3 i  terms  which  appear  in  the  current  variable  formulation. 

However,  available  data  and  existing  models  for  diodes,  transistors  and 
SCR's  are  given  in  terms  of  current  and  voltage  variables  instead  of  flux 
linkages  and  charge  variables.  This  and  the  need  to  use  a  program  such  as 
SPICE  or  SCEPTRE  for  the  electronic  component  modeling  dictated  the  use  of 
current  variables  for  modeling  the  generator.  SCEPTRE  normally  requires 
that  currents  be  the  state  variables  in  inductor  elements. 

2.4  FORMULATION  FOR  NON-LINEAR  INDUCTANCE 

Saturation  effects  cause  the  inductance  coefficients  to  be  non-linear 
functions  of  the  machine  currents.  The  circuit  model  equations  for  direct- 
phase  current  formulation  with  saturation  effects  are  now  developed.  The 
terms  arising  due  to  dependence  of  inductance  on  the  currents  are  illustrated. 

The  generator  circuit  model  equations  for  direct-phase  current  formula¬ 
tion  and  inductance  as  a  function  of  current  are 

V  =  ri+l3T*<3T><i>  <6> 
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In  the  second  term  on  the  right-hand  side  of  the  equation  the  induc¬ 
tances  are  those  defined  in  the  appendix  by  Equation  (50)  except  that  the 
L  and  M  coefficients  are  modified  by  saturation.  The  third  term  involving 
dL/dt  will  be  treated  subsequently. 

The  assumptions  made  on  the  effect  of  saturation  on  the  inductance  terms 
are  discussed  in  Appendix  A.  Equation  (53),  repeated  here,  illustrates  how 
one  of  the  inductance  terms  is  a'ected  by  saturation. 


'aas 


=  C  Le  +  C,.L_  cos  20 
aa  s  aa  m 


(53) 


^aa*"aa 

Implied  here  is  that  the  shape  of  the  inductance  curve  as  a  function  of 
rotor  position  is  not  changed  by  saturation,  (the  validity  of  this  assump¬ 
tion  is  discussed  in  Appendix  A).  This  form  further  implies  that  saturated 
inductance  measurements  can  be  made  on  the  direct  axis  as  a  function  of  net 
excitation  contributed  by  the  various  currents.  The  net  excitation  current 
in  terms  of  equivalent  main  field  current  is 

'x  '  *F  t  Vo  +  NFa  t'a  C0S(B)  +  'b  cos(e  '  X1 

♦  ic  cos (a  +  ^)]  (7) 


To  determine  the  C  ,  the  rotor  is  aligned  with  the  "a"  phase  axis  and 

ad 

Laa  is  measured  as  i  is  varied  over  an  appropriate  range  giving  a  curve  as 
shown  in  Figure  1. 

The  equation  for  Lga  as  a  function  of  i  is  obtained  by  appropriate 
polynomial  curve  matching  techniques  and  results  in  the  following  form. 


‘aas 


*  Lao<>  +  *1  fx  *  Vx  +  Vx  + 


(8) 


LaoCaa 


“aas 


Thus  C 

aa 

with  1  . 


is  the  normalized  polynomial  of  the  variation  of  the  inductance 
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I  FIGURE  1.  VARIATION  OF  L  ,  WITH  1  (SATURATION  EFFECTS) 

da  X 


C  terms  are  required  for  each  of  the  equations  listed  in  Equation  (50) 
Measurements  are  needed  to  determine  the  following  C  polynomials:  Caa,  Cflp 
CaDt  C3q,  Cpp,  Cp0,  CDD.  These  requirements  are  listed  in  more  detail  in 
a  later  section. 

These  saturation  coefficients  are  associated  with  the  inductance  terms 
as  shown  below  in  Equation  (9)  where  the  s  subscript  is  dropped 


La,  ■  Caa  <Ls  +  Cn,  c°s<2e» 

L.b  *  C.a  <-Hs  +  Lm  C°s<2e  '  T»  '  Lb» 
Lac  -  caa  K  *  Lm  Cos(2e  *  T»  *  Lc. 
LaF  s  CaF  *MaF  Cos^0^  *  LFa 


L-n  -  C_n  (M,n  Cos(q)) 


LaQ  "  CaQ^MaQ  sin(®))  “  lq8 
Lbb  ■  Ca.<Ls  *  U  c“<29  +  T» 

Lbc  ■  Ca.<-Ms  +  L»  Cos<26»  '  Lcb 
LbF  *  <v<V  Cos<9  *  T»  *  LFb 
LbO  *  CaD^MaD  c°*(e  -  T»  ■  LDb 
LbQ  *  caq(MaQ  si"(e  -  T»  "  «-Qb 
Lcc  '  caa(Ls  +  L«  Cos<29  -  X» 

LcF  =  caF(HaFCos‘6+T))  *  LFc 
lcD  ’  caD<MaD  c°s<e  +  T»  ■  L0c 
UcQ  *  CaQ(MaQ  s1n(6  +  T»  *  LQc 


lff  = 

cfflf 

lfd  = 

cfdmr 

3  ldf 

lfq  = 

0 

=  lqf 

ldd  3 

cddld 

ldq  3 

0 

3  lqd 

lqq  =  lq 

Next  consider  the  third  term  on  the  right  hand  side  of  Equation  (6), 
(dL/dt)(I). 
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Consider,  for  example,  the  first  row  of  the  matrix  product  (dL/dt)(I). 


Laa><a 


+ 


Lab}ib  +  Lac)ic  +  LaF)iF 


+  (^‘LaD)iD  +  (^’LaQ)iQ 


(10) 


now  noting  that  the  inductances  are  functions  of  the  net  excitation,  ix,  and 
i  is  given  by  Equation  (7),  the  derivative  in  the  first  term  of  Equation  (10) 
expanded  becomes. 


a  '  3L,,  3iv  3L,_  .... 

d  ■  _  3d  X  ^  33  39 

dt  Laa  ”  3ix  3t  30  3t 


(11) 


3L 


aa 


3i. 


3ip 

3t 


31, 


+  N 


FD  3t 


+  NFa'irf)Ct,s<e>  *  HFaSa  si"<8> 


36 

3T 


+  NFa<TT>C°s<9 


T>  -  "Fa'b  sin(9  ‘  T~>  H 


*  NFa<Tf>Cos<9  +  T>  -  NFa  S1n<e  +  T>  If 


3L 


aa  36 
36  dt 


(12) 


The  expressions  for  the  other  terms  after  expanding  the  derivative  terms 
of  Equation  (10): 

dt  Lab’  3t  Lac’  dt  LaF’  dt  LaD’  dt  uaQ’ 

will  be  the  same  as  Equation  (12)  except  for  the  quantities  outside  the 
square  brackets  where  the  inductance  in  question  will  appear.  When  these 
terms  are  all  expanded  and  terms  are  collected,  the  first  row  of  the 
L  jj-  +  (^)  I  matrix  becomes 
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■  CVf,  c°s<9>  +  Lscaa  +  Lmcaa  c«Ue)3  -jf 

*  [AaNfa  Costs  -  %■)  -  MsCm  ♦  LmCaa  Cos (26  -  ^)]  ^ 

*  tAaHFa  Cos*e  +  t'  ■  Mscaa  +  LmCaa  Cos<2e  +  t*3  ~3T 

*  tAa  +  CaFMaF  Cos  (9^  IT 

+  tAaNFD  *  caDMaD  Cos  t 0 )  3 
+  [CaQMaQ 

-  [AaB  +  2caaCaLm+(CaFMaFiF  +  CaDiD)si"(e)-CaqV'Q  Cos«8>3  at 


B  *  Npa[(a  Sin(6)  *  ib  S1n(«  -  f-)  *  (c  Sin{8  +  f)] 


Ca  =  ia  Sin(2e}  +  ib  Sin(2e  -  +  1C  Sin(2e  +  Q) 

Similar  expansions  are  required  for  each  of  the  rows  of  the  matrix  terms 
corresponding  to  the  terms 


3  A.  3A_  3Ar  3A 


OAb  OAc  OAF  °AD  A  "*0 

at  •  at  »  at  •  at  »  and  TF 


The  complete  listing  of  these  terms  is  included  in  Appendix  B. 

Further,  one  must  note  and  include  the  saturation  effects  on  the 
terms.  Again  using  phase  "a"  terms  as  example 

Laa  ‘  Caa<Ls  *  Lm  C°s<29» 


and  noting  from  Equation  (8) 


aa 


a  1  +  Vx  +  *2'l  +  *£  + 


(15) 


*  (L$  +  Lm  Cos(2e) ) (a-j  +  2a2ix  +  3a31*  +  •••)  (16) 

Again  the  derivative  of  the  other  inductance  terms  with  respect  to  ix 
may  be  written  in  a  similar  fashion. 

Finally  the  derivatives  of  the  inductance  terms  with  respect  to  e  are 
obtained.  The  "a"  phase  terms  are  listed 

TT  *  -  2  caaLn,  S,"<2|» 

if  '  -  2  Cao*'ra  Si"<29  -  T> 

if  *  -  2  CaaLin  s'"<2e  *  T> 

~TT  *  "  caF^aF  s1"*6' 

TT  *  -  CaDMaD 

1?  *  •  caQMaQ  Cos(6>  •  (,7> 

2.5  SCEPTRE  EQUIVALENT  CIRCUIT 

The  circuit  model  equations  developed  in  previous  sections  are  now  in  a 
form  from  which  the  equivalent  circuit  may  be  readily  discerned.  An  equivalent 
phase  "a"  circuit  that  is  easily  analyzed  by  the  SEPTRE  program  is  shown  in 
Figure  2.  Equations  for  the  other  circuits  of  Figure  20  are  developed  in  a 
similar  manner  and  equivalent  circuits  built  for  them.  A  program  in  SCEPTRE 
can  then  be  written.  Such  a  program  is  listed  in  Appendix  C. 
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in 


si 


71 


61 


51 


41 


31 


21 


1 


'  tVro  4  c.d\o  «»*«•)]  -jt 

dtc  autual 

-  ♦  C.J  f  Cos(»)J  -a*-  >  Inductance 

*  "  term 

-  [»,Nf,  Cos(4  ♦  y*  ♦  lwl  -jf 
*  [A»nF»  Co,(*  '  T*  4  W  IT 

-  c.oVq  Cos(,)  & 

^.o'o  s'nt*>  ar 


C4F,Wf  If 
*..yc  S1"<2*  4  T->  3F 
*..Vb  *'"*•  -  T>  M 

2£„L,  SIn(2»)  |f 


1\"f.  <*,(.)  ♦  L,„] 


FIGURE  2.  EQUIVALENT  CIRCUIT  FOR  PHASE  *•*  OF  SYNCHRONOUS  MACHINE. 


2.6  INCLUSION  OF  PRIME-MOVER  AND  EXCITATION  CONTROL  MODEL 


While  SCEPTRE  was  primarily  developed  to  facilitate  network  analysis, 
it  can  be  used  to  analyse  other  dynamic  systems.  Particularly,  it  has  the 
ability  to  analyze  systems  for  which  transfer  functions  are  available. 

SCEPTRE  requires  that  the  transfer  function  be  converted  to  state  equations 
and  the  state  equations  entered  into  the  program.  (The  details  of  this 
conversion  procedure  are  given  in  the  manual  (Reference  2).)  This  capability 
of  SCEPTRE  can  be  used  to  advantage  in  simulating  the  prime-mover  charac¬ 
teristics  of  the  system.  The  prime-mover  model  is  converted  to  a  suitable 
program  and  the  output  quantity  of  interest  is  the  mechanical  torque,  Tm  . 
quantity  becomes  an  input  for  the  synchronous  machine  model  through  Equation 
(3),  which  is  the  equation  for  rotor  angle  acceleration. 

Similarly,  excitation  control  system  models  may  be  added  to  the  machine 
simulation.  Excitation  control  is  generally  based  on  monitoring  the  terminal 
voltage  of  the  generator  and  changing  the  excitation  field  voltage  in  some 
fashion  as  the  response.  The  control  system,  modeled  either  as  a  transfer 
function  or  an  electronic  network,  is  easily  added  to  the  present  machine 
model,  since  in  the  program  both  the  terminal  voltage  and  the  excitation 
voltage  are  accessible  variables.  The  present  work  has  not  been  concerned 
with  modeling  either  the  prime-mover  or  the  excitation  voltage  control 
scheme,  and  the  above  observations  are  offered  only  to  show  how  these  models 
can  be  incorporated  in  the  present  simulation. 
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SECTION  III 


THE  TRANSFORMER  MODEL 


3.1  OBJECTIVE 

In  this  section  various  problems  and  options  are  considered  which  will 
determine  the  mathematical  and  computer  program  of  a  three  phase-three  leg 
core  transformer  with  wye-wye  or  delta-wye  connections.  The  model  is  to 
be  adequate  for  handling  unbalanced  transient  loads,  start  up  transients 
for  the  system,  and  include  saturation  conditions. 

The  first  year's  work  on  the  transformer  model  has  involved  the  trans¬ 
former  equivalent  circuit,  the  magnetic  core  characteristics,  relation  of 
circuit  parameters  to  magnetic  core  characteristics  and  study  of  phenomena 
that  may  be  potential  problems  such  as  initial  inrush  of  current  and  ferro- 
resonance. 


3.2  THE  TRANSFORMER  EQUIVALENT  CIRCUIT 


The  equation  defining  the  circuit  relationship  for  each  transformer 
winding  is: 


where 


dx 


vkj  =  R 


kj  nkj  +  dt 


J£ 


k 

j 

xkj 


a,  b,  c  for  the  three  phases 

1 ,  2  for  primary  and  secondary 

is  the  terminal  voltage  of  the  kj  winding 

is  the  current  in  the  kj  winding 

the  linkages  of  winding  kj  . 


(18) 


The  interrelationship  of  the  currents  and  voltages  depend  on  the  network 
connections.  These  relationships  are  automatically  imposed  by  the  SCEPTRE 
modelling  program. 

^*ki 

The  relates  the  electric  network  equations  to  the  magnetic  field 
of  the  three  phase-three  leg  core  of  the  transformer. 


In  general  the  linkages  of  any  winding  are  some  non-linear  functions  of 
all  of  the  winding  currents  as  in  Equation  (19) 

Aki  =  xkj^al  *  V  ’  1bl  ’  ^b2  *  jcl  *  W* 

Due  to  hysteresis  effects.  Equation  (19)  is  not  a  single  valued  function 
but  the  value  for  x^,.  depends  not  only  on  the  particular  state  of  the  six 
current  variables  but  also  on  whether  each  of  the  six  is  increasing  or 
decreasing.  A  precise  representation  of  a  three  phase  -  three  leg  core 
transformer  by  Equation  (19)  would  require  an  unreasonable  amount  of  data 
whether  obtained  by  laboratory  measurements  or  by  calculations  using 
magnetic  core  data. 

Subsequent  paragraphs  in  this  section  discuss  some  assumptions  commonly 
made  which  reduce  the  amount  of  data  required  to  model  the  magnetic  charac¬ 
teristics.  Some  of  these  approximations  or  assumptions  leave  one  with  a 
reasonably  accurate  model  because  of  conventional  construction  practices, 
others  represent  a  trade-off  between  the  amount  of  data  required  and  the 
accuracy  of  the  model. 

The  results  of  assuming  a  constant  leakage  inductance  and  constant 
turns  ratio  between  windings  of  the  same  phase  are  illustrated  in  the 
development  of  Equation  (20)  through  Equation  (29). 

Define  mutual  and  leakage  flux  in  phase  "a"  by  Equations  (20)  and  (21) 


similar  expressions  apply  to  the  ”b"  and  "c"  phase  windings. 

In  terms  of  linkages  with  an  effective  primary  turns,  Nj  Equation  (18) 
can  be  rewritten  for  phase  "a"  primary  and  secondary 

'a,  ■  Ral  <«1  ♦  L»al  TT  +  IT 


V.-  =  i.-,  +  sar-  + 


N7  dxam 
c  am 


a2  a2 


ta2  sr T  Ny-  -jr 
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If  Xan)  is  understood  to  be  in  terms  of  primary  turns.  Thus  in  a 
three  phase  transformer  there  are  three  linkage  variables  xam ,  x^  and  xcm  . 

Further  implied  by  the  constant  turns  ratio  of  a  given  phase  is  that 
there  are  three  excitation  variables  instead  of  six. 


Fa  =  Nai  ^1  +  Na2  ^2 
FB  =  Nbl  1bl  +  Nb2  ^2 
Fc  =  Ncl  icl  +  Nc2  ic2 

Thus  one  has  the  linkages  as  functions  of  three  excitations 


(24) 

(25) 

(26) 


xkm=  WFa  ’  Fb  ’  Fc^ 

and  the  derivative  term  of  Equation  (22)  becomes 


(27) 


dXam  3*,m  dFa 

am  =  ani  _ _ a 

dt  aF_ 


axam  dFk 
.  am  b 

dt  3Fb 


ax,m  dF 

..  ,  am  c 

HT  TF^T  dt 


(28) 


Thus  the  complete  three  phase  magnetic  field  to  circuit  relationship 
can  be  expressed  in  matrix  form 


dx,„ 

3X,,, 

ax. 

ax 

dF, 

am 

am 

am 

am 

a 

dt 

3Fa 

3Fb 

3Fc 

"3T 

dlbn, 

3xbm 

3Xbm 

3xbm 

dFb 

dt 

af 

aF. 

nr 

a 

b 

c 

dx 

ax„„ 

ax 

ax 

dF. 

cm 

cm 

cm 

cm 

c 

dt 

[3Fa 

3Fb 

3FcJ 

nr 

(29) 


Equation  (29)  implies  that  nine  magnetic  characteristics  are  required. 
However,  each  of  these  are  families  of  hysteresis  loops. 

Note  that  each  term  in  the  square  matrix  applies  to  both  primary  and 
secondary  of  the  appropriate  phase.  Further,  if  constant  leakage  inductance 
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had  not  been  assumed  this  matrix  would  be  a  6  x  6  requiring  36  magnetic 
characteristics.  Otherwise  the  analysis  would  be  the  same. 


3.3  THE  MAGNETIC  CORE  DATA 

The  major  tasks  in  applying  the  circuit  equations  (Equation  22),  are 
in  obtaining  the  data  required  for  the  incremental  inductances  of  Equation 
(29),  storing  the  large  amount  of  data  in  the  computer  in  an  efficient  manner 
and  developing  the  logic  to  access  the  data  at  the  appropriate  operating 
point. 

The  type  and  amount  of  data  are  discussed  briefly  in  Nakra  and  Barton 
(Reference  10)  for  a  somewhat  simpler  system.  As  a  method  of  storing  B-H 
curve  data,  several  mathematical  forms  have  been  used.  Macfadyen  et  al 
(Reference  13)  discusses  representation  of  magnetization  curves  by  exponen¬ 
tial  series.  Others  use  hyperbolic  forms  or  polynomials.  Manly  (Reference  9) 
compares  several  methods.  Somewhat  less  accurate  but  more  efficient  is 
the  method  given  in  the  Air  Force  report  AFAPL-TR-76-102  (Reference  14). 

In  terms  of  the  B-H  curve  instead  of  the  X-i  curve  the  equations  are 

B  =  Bs(H  +  KHc)/[Hc(Bs/Br  -  1 )  +  |H  +  KHj]  -  KBQ  (30) 

where 

Bs  is  the  saturation  value  of  flux  density  for  the 
major  loop 

Br  is  the  residual  value  of  flux  density  for  the 
major  loop 

Hc  is  the  coercive  force  for  the  major  loop 

K  is  +  1  to  represent  the  lower  and  upper  curves. 

The  Bq  allows  representation  of  the  minor  loops  (and  hence  the  uppers 
and  downers)  where  the  minor  loops  intersect  at  point  (+  Hm  ,  +  Bm)  and  BQ 
is  given  by 


Bc(H+HJ 

s '  m  c 


BfH  -Hj 
s  m  c 


Rc<VBr-,>t|Hm*Hcr+Hc(Bs/Br-l)+|Hin-Hc| 


(31) 


Figure  3  snows  curves  computed  by  these  formulas. 
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Thus  only  data  at  several  points  on  the  major  loops  are  required,  the 
Br  ,  Bs  and  Hc  values. 

The  B  and  H  variables  are  cornnon  design  variables  but  the  X  and  1 
variables  are  more  common  in  measurements  and  are  needed  in  network  equa¬ 
tions. 


3.4  A  SIMPLIFYING  APPROXIMATION 


tach  entry  in  the  incremental  inductance  matrix  of  Equation  (29)  re¬ 
quires  the  identification  of  an  operating  point  in  terms  of  the  net  excita¬ 
tion  of  the  combined  Fa  ,  F^  and  F£  .  For  each  pair  of  values  for  Fjj  and 
Fc  *  Fa  varied  over  its  range  to  determine  a  X^  major  or  minor  hysteresis 
loop.  (For  each  incrementing  of  either  F^  or  Fc  requires  a  new  family  of 
hysteresis  loops.)  By  using  Equations  (30)  and  (31),  data  for  major  loops 
only  are  required.  This  still  is  a  large  amount  of  data  to  collect,  store 
and  access. 


A  slightly  different  approach  will  allow  some  approximations  which 
will  greatly  reduce  the  amount  of  data  required  and  make  definition  of  the 
operating  point  easier. 

Returning  to  Equation  (27)  and  (28)  let  us  define  a  net  excitation  of 
phase  "a" 


F;  =  Fa 


+  Wb 


+  TacFc 


then 


(32; 


dFa  dFh  dFr 

■JT  +  Tab  HT  *  Tac  ~ST 


(33) 


Equation  (33)  requires  one  incremental  inductance  term  and  one  set  of 
hysteresis  loops  xam  vs  Fa  as  shown  in  Figure  4.  However,  the  apparent  turns 
ratios  Tat)  and  Tac  will  be  variables  depending  on  specific  levels  of  Ffl  ,  F^ 
and  Fc  . 

Now  if  one  can  neglect  any  change  in  shape  of  the  xflm-  Fj  curves  for 
differing  levels  of  F^  and  Fc  ,  then  F^  and  F£  are  accounted  for  by  horizontal 
shifts  in  the  loops  of  Figure  4.  At  each  point,  the  shift  of  the  vertical 
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axis  from  the  plot  of  xam  as  a  function  of  F#  alone,  gives  a  measure  of  T#b 
and  TflC .  (A  set  of  major  loops  varying  F#  over  the  full  range  for  the  set 
of  conditions  defined  in  Table  1  can  be  used  with  multiple  regression  to 
define  rab(Fa  .  Fb  ,  Fc)  and  Tac(Fa  ,  F„  ,  F,.).) 

Thus  one  family  of  loops  and  equations  for  Tflb  and  Tac  would  be  all 
the  data  required  for  phase  "a”  magnetic  characteristics.  Similar  data  and 
representation  would  be  required  by  phases  "b"  and  "c”  windings. 

3.5  DETERMINING  THE  OPERATING  POINT 

In  a  computer  program  which  performs  a  numerical  integration  to  de¬ 
termine  the  system  transient  solution,  a  subroutine  is  required  to  de¬ 
termine  the  incremental  inductance  at  each  operating  point,  ihis  incre¬ 
mental  inductance  is  dependent  on  the  state  of  each  of  the  current  variables 
and  further  the  incremental  inductance  is  dr pendent  on  whether  each  variable 
is  increasing  or  decreasing.  The  use  of  excitation  defined  by  Equation  (32) 
and  interpolation  between  curves  such  as  those  shown  in  Figure  4  would 
determine  an  approximate  incremental  inductance. 
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TABLE  1 


SECTION  IV 

THE  AC  RESONANT  CHARGING  CIRCUIT 


4.1  OBJECTIVE 

This  part  of  the  project  was  the  development  of  the  ac  resonant  charging 
circuit  model  which  is  to  be  the  load  on  the  generator  and  transformer  dis¬ 
cussed  in  earlier  sections.  The  characteristics  of  the  series  RLC  resonant 
charging  circuit  are  well  known  (Reference  22)  and  the  sequence  of  switching 
events  defined  by  prior  knowledge  of  the  particular  application.  The  work 
on  this  part  of  the  project  thus  consists  of:  selection  of  an  SCR  model  that 
will  adequately  represent  turn-on  and  turn-off  characteristics  at  the  power 
level  and  with  the  resonating  inductors  of  this  particular  application; 
assure  that  the  model  works  when  using  series  and/or  parallel  connected 
SCR's  for  higher  power  applications;  convert  the  model  to  a  form  compatible 
with  the  integrated  system  model. 

The  accomplishments  reported  here  are  the  selection  of  the  SCR  model 
with  appropriate  modifications  for  this  application,  testing  of  the  model 
using  CAD  program  SPICE2,the  conversion  of  a  SPICE  2  SCR  model  to  a 
SCEPTRE  model,  and  initial  afforts  of  ac  resonant  charging  using  SCEPTRE. 

4.2  SCR  MODEL  SELECTION 

Five  options  were  developed  and  considered  as  possible  approaches  to 
the  objective  of  the  SCR  model  selection.  A  brief  summary  of  each  option 
is  as  follows: 

a.  Develop  a  Model  in  House 

This  option  was  viewed  as  not  well  suited  to  the  time  frame  of  the 
project.  It  would  be  used  only  in  the  event  that  no  suitable  existing 
model  could  be  found. 

b.  The  Two  Bipolar  Junction  Transistor  (2-BJT)  Model  (Reference  23) 

This  model  became  the  favored  candidate  because  of  the  following.  It 
was  a  long  established  model  with  a  great  deal  of  available  supporting  litera¬ 
ture.  The  model  had  been  verified  by  University  of  Callfornia-Berkeley 
(Reference  23)  to  adequately  simulate  the  switching  and  operational 
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characteristics  of  an  SCR.  A  method  for  determining  model  parameters  from 
manufacturers  specification  sheet  data  had  also  been  developed  by  University 
of  California-Berkeley.  This  model  had  already  been  demonstrated  using 
SPICE  2  and  appeared  readily  adaptable  to  other  advanced  CAD  programs  such 
as  SCEPTRE  (Reference  23). 

c.  The  Basic  Three  Junction  (Intrinsic)  Model 

The  intrinsic  SCR  model  is  fundamentally  similar  to  the  2-BJT  Model. 

No  advantages  over  the  2-BJT  Model  were  readily  apparent.  One  disadvantage 
was  apparent.  That  being  an  increased  difficulty  to  the  operator  for  use 
with  SPICE  2. 

d.  Expanded  Three  Junction  (J.  Bowers)  Model  (References  24,30) 

This  model,  developed  by  J.  Bowers  at  the  University  of  South  Florida, 
is  a  complex  model  which  is  an  expandion  of  the  basic  3-Junction  Model 
which  again  is  fundamentally  similar  to  the  2-BJT  Model.  J.  Bowers  Model 
was  viewed  as  too  complex  for  the  application  desired.  In  addition  the 
complete  model  parameters  are  very  difficult  to  obtain.  On  the  other  hand, 
this  model  could  be  simplified  to  more  nearly  resemble  the  Basic  3-Junction 
Model.  The  University  of  California-Berkeley  method  of  parameter  determina¬ 
tion  could  then  possibly  be  extended  to  this  model.  This  model  has  also 
been  extensively  verified  in  performance  using  SCEPTRE.  A  difficulty  was 
apparent,  however,  in  the  case  of  adaptation  for  use  with  another  CAD 
program  such  as  SPICE  2. 

e.  Georgia  Tech  Terminal  Characteristics  Model  (Reference  25) 

This  option  offered  the  possibility  of  a  high  power  SCR  model  as  called 
for  by  the  research  project  at  hand.  Unfortunately,  development  of  the 
model  was  only  partially  complete  and  therefore  could  not  be  fully  evaluated. 
To  the  extent  that  development  had  proceeded,  it  appeared  that  determination 
of  model  parameters  might  be  difficult. 

There  were  several  important  considerations  in  selecting  the  model. 

Among  these  are  first,  the  ease  of  obtaining  numerical  values  for  model 
component  parameters  from  manufacturer's  data  or  by  measurement,  second, 
the  integrity  of  the  model  as  established  by  the  literature  or  by  test.  The 
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model  must  adequately  represent  features  important  in  this  particular  applica¬ 
tion.  Third,  in  this  project  the  compatability  of  the  model  with  CAD  programs 
such  as  SPICE  2  and  SEPTRE. 

The  2-BJT  model  shown  in  Figure  5  was  altered  slightly  and  used  in 
this  project. 

Appendix  E  shows  the  development  of  the  BJT  as  it  evolved  from  the 
Ebers-Moll  model  through  the  work  of  Ki  (Reference  23)  and  Getreu  (Reference 32). 
Appendix  E  shows  the  development  of  the  SCR  model  from  the  BJT,  the  determina¬ 
tion  of  the  element  values  required  for  SPICE  2  modeling  of  the  SCR  and  the 
determination  of  the  element  values  required  for  the  SCEPTRE  modeling  of  the 
SCR.  The  SCEPTRE  values  are  obtained  as  conversions  from  the  more  familiar 
SPICE  2  quantities. 

4.3  SCR  MODEL  IN  AC  RESONANT  CHARGING 

The  CAD  program  SCEPTRE  is  to  be  used  in  modeling  the  composite  system 
of  the  ac  generator,  transformer,  and  the  load  consisting  of  ac  or  dc 
charging  circuit.  While  it  is  recognized  that  SCEPTRE  is  generally  slower 
than  SPICE  2,  SPICE  2  does  not  have  the  capabilities  required  for  modeling 
the  generator  non-linear  time  varying  parameters.  However,  SPICE  2  is  faster 
and  more  convenient  to  use  in  electronic  circuits  such  as  the  ac  resonant 
charging  circuit.  In  the  following  work  the  model  development  and  testing 
of  the  SCR  in  ac  resonant  charging  is  done  using  SPICE  2,  with  the  conversion 
to  SCEPTRE  format  as  a  later  task. 

a.  SPICE  2  Single  Loop  Analysis 

A  SPICE  2  program  was  written  and  ran  for  the  circuit  of  Figure  6.  The 
results  of  the  SCR  ac  resonant  charging  circuit  single  loop  analysis  revealed 
an  oscillation  following  commutation  (Ref.  to  Fig.  7).  The  oscillation 
appears  to  result  from  the  non-linear  resonant  behavior  of  the  junction 
capacitances  with  the  circuit  inductance.  This  is  supported  by  the  fact 
that  removal  of  the  inductor  from  the  circuit  resulted  in  normal  commutation 
with  no  oscillation  after  SCR  cut-off.  The  oscillation  may  be  suppressed  by 
addition  of  a  snubber  to  the  circuit  (Ref.  to  Fig.  8);  however,  investiga¬ 
tion  of  the  fundamental  nature  of  the  oscillation  is  continuing. 


FIGURE  6.  SPICE  CIRCUIT  DIAGRAM  FOR  THE  SINGLE  LOOP  AC  RESONANT  CHARGING  SYSTEM 


FIGURE  7.  SPICE  2  SIMULATION  OF  THE  AC  RESONANT  CHARGING  CIRCUIT  OF  FIGURE  6. 
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FIGURE  7a.  SPICE  2  PROGRAM  LISTING  FOR  CIRCUIT  OF  FIGURE  6  WITH  SNUBBERS 
SIMULATION  RESULTS  ARE  SHOWN  IN  FIGURE  4.5. 


SIMULATION  OF  SINGLE  LOOP  AC  RESONANT  CIRCUIT  WITH  A  SNUBBER 
9).  THE  "MODIFIED  HU-KI  MODEL"  WAS  USED. 
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FIGURE  9.  SINGLE  LOOP  CIRCUIT  WITH  SNUBBER. 


b.  SPICE  2  Two  Loop  Analysis 


A  second  SCR  and  charging  capacitor  loop  was  added  to  the  circuit  as 
shown  in  Figure  10.  The  associated  SPICE  2  Input  listing  is  shown  in 
Figure  11. 

The  two  loop  analysis  revealed  a  problem  with  the  second  branch  SCR 
turning  on  prior  to  application  of  the  gate  2  trigger  signal.  This  early 
firing  of  the  second  branch  SCR  was  suspected  to  be  related  to  the  oscilla¬ 
tion  experienced  following  turn-off  of  the  first  branch.  Accordingly,  a 
program  was  run  in  which  all  things  were  kept  the  same  with  the  exception 
that  no  gate^  pulse  was  applied  to  SCR^  .  The  circuit  operated  correctly 
in  that  SCRg  did  not  turn-on  until  application  of  the  gateg  signal.  Snubber 
circuits  were  then  added  in  an  effort  to  and  suppress  the  early  turn-on. 

This  effort  has  not  yet  been  successful. 

In  summary,  the  efforts  were  inconclusive  and  effort  is  still  under 
way  to  solve  problems  incurred.  The  problems  are  that  the  circuit  inductance 
has  a  resonant  interaction  with  a  SCR  junction  capacitance  which  causes 
unacceptable  circuit  oscillation  following  turn-off  of  the  device.  The 
circuit  oscillation  then  causes  uncontrolled  turn-on  of  other  SCRs  in  the 
circuit.  Conventional  snubber  circuits  have  not  prc/en  successful  in 
suppressing  this  problem. 
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FIGURE  II.  SPICE  2  TNO  LOOP  AC  RESONANT  CHARGING  SYSTEM  INPUT  DECK 
FOR  CIRCUIT  SHOWN  IN  FIGURE  10. 


SECTION  V 

CHOICE  OF  VARIABLES 


5.1  OBJECTIVE 

When  solving  systems  of  equations  using  numerical  methods  and  in  par¬ 
ticular  when  the  equations  are  non-linear,  the  choice  of  system  variables 
may  have  significant  effect  on  the  computation  time  and  on  the  reliability 
of  the  results.  In  this  section  the  nature  of  these  effects  are  examined 
even  though  the  decision  to  use  CAD  programs  such  as  SPICE  2  and  SCEPTRE 
may  preclude  such  considerations. 

5.2  CHOICE  OF  VARIABLES 

Chua  and  Lin  (Reference  11)  have  shown  that  the  selection  of  a  certain 
set  of  state  variables  depends  on  the  nature  of  the  nonlinear  elements.  For 
a  nonlinear  capacitor  characterized  by  a  nonmonotonic  voltage-controlled  q-v 
curve,  for  example,  the  capacitor  voltage  must  be  chosen  as  the  state- 
variable.  Similary,  inductor  current  must  be  chosen  as  the  state  variable 
for  a  nonlinear  inductor  characterized  by  a  nonmonotonic  current-controlled 
A-i  curve.  However,  if  the  capacitor  and  inductor  characteristic  curves  are 
strictly  monotonic,  then  either  capacitor  voltage  or  charge  may  be  chosen  as 
the  capacitor  state  variable,  and  for  the  inductor,  the  state  variable  may 
be  choseii  «;s  either  inductor  flux  linkage  or  current. 

However,  it  is  shown  (Reference  12)  that  it  is  advantageous  to  choose 
capacitor  charge  and  inductor  flux  linkage  as  the  state  variables  when  a 
numerical  integration  algorithm  is  used  because  this  particular  choice  of 
state  variables  reduces  the  global  error.  Global  error  is  the  error  accrued 
over  the  finite  time  interval  during  which  the  numerical  integration  is 
performed. 

To  see  why  this  is  so,  consider  a  nonlinear  system  (Reference  12) 
modeled  by  the  equations 

x-  f(x ,  t)  x(0)  =  (34) 

where  the  bar  implies  the  use  of  vectors  or  matrices. 
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To  estimate  the  growth  of  errors  over  some  time  interval ,  say  time  t  *  0 
to  t  =  tf  when  solving  Equation  (34)  numerically,  consider  the  differential 
equation  which  results  when  the  initial  condition  is  perturbed  so  that 
£.(0)  £©  +  5*(0)  =  jSq  +  '  ™en  ^rom  ^uatlon  (34)  the  perturbed  equa¬ 

tions  become 

(x(t)  +  6xU))  ■  f(x(t)  +  6x(t),  t)  (35) 

6x(0)  =  ^ 

Assuming  6x(t)  is  small  and  expanding  the  right  side  of  Equation  (35) 
by  a  Taylor  series  gives 

^  (x(t)  +  6x(t))  s  f(x),  t)  +  |=  6x(t)  (36) 

6x(0)  *  ^ 

Subtracting  Equation  (36)  from  Equation  (34)  gives 


.  ,  3f 

3x  *  (~J  6x 


6x(0) 


where  the  matrix  (3f/3)0  as  a  time-varying  Jacobian  matrix  Jx(t).  Equation 
(37)  represents  a  linear  time-varying  system;  thus  no  simple  solution  exists. 
In  the  case  where  x  and  f  are  scalar  quantities,  the  solution  to  Equation  (37) 
is 


6x(t)  =  «0  e 


/  dJOdT 
o  * 


Equation  (38)  indicates  that  whenever  Jx(t)  >  0  ,  the  error  (accumulated 
at  time  t  in  the  exact  solution  to  the  scalar  version  of  Equation  (34),  due 
to  an  error  in  the  initial  condition)  actually  decreases  exponentially 
with  time.  Since  the  concern  is  the  accrued  effects  of  errors  made  at  all 
time  steps  (and  not  just  that  due  to  an  error  in  the  initial  condition), 
proceed  to  determine  the  accrued  error  E  *tf,  ie.,  the  sum  of  errors  from 


time  t  =  0  to  time  t  =  t4 


To  this  end  it  is  assumed  that  the  error  per 
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unit  time  (which  Is  denoted  by  c^.)  Is  constant.  Then  the  error  made  In  an 
Interval  of  time  dt.  Is  dt^  .  From  Equation  (38),  the  growth  of  error 
for  t-j  <  t  Is  given  by 

t 

/  JytOdT 

tl.t  t,  X 

e  =  dt'j  e  (39) 


0»tf 

so  that  the  accrued  error  E  at  time  t^  is  the  sum  of  the  error  originating 
for  all  times  t^  <  t^  ,  so  that 


(40) 


To  demonstrate  error  growth,  Calahan  (Reference  11)  used  the  circuit  of 
Figure  12,  where  the  nonlinear  capacitance  is  strictly  monotonic.  If  the 
state  equation  is  written  in  terras  of  the  capacitor  voltage  we  have 

<CT+Cd  e*V)  ar=  -  if  •  isiex',-D;  ‘  -  «  m 

The  Jacobian  is  given  by 


where 


JV  = 


ACd  eXV(Vi)  +  Is(eXV‘1) 
(CT  +  C°  eXv)Z 


(ff  +  AISCXV) 

(CT  +  C°  eXv) 


A  =  40  . 


(42) 


During  switching,  when  vs  <  0  and  v  >  0,  the  first  term  shifts  the 
eigenvalue  Jy  toward  the  right  half  plane.  Indeed  if  v$  is  made  suffi¬ 
ciently  negative,  the  eigenvalue  Jy  can  be  made  as  positive  as  desired. 
Thus,  during  this  time  when  Jy  is  positive,  truncation  errors  introduced  by 
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the  integration  formula  used  will  be  amplified  according  to  Equation  (40). 

By  choosing  charge  as  the  state  variable,  the  resulting  eigenvalue  Jp  can 
be  shown  to  lie  in  the  left  half  plane  during  switching,  so  that  truncation 
errors  will  decay  according  to  Equation  (40).  Finally,  it  should  be  noted 
that  in  the  case  of  linear  capacitors  or  inductors,  either  q  or  v  for  a 
capacitor,  or  X  or  i  for  an  inductor  may  be  chosen  as  the  state  variable 
(Reference  11). 

5.3  SCEPTRE  VARIABLES 

In  Section  VII  the  decision  to  use  the  CAD  system  called  SCEPTRE  is  dis¬ 
cussed.  Whether  SCEPTRE  or  the  other  system  considered  SPICE  2,  is  used, 
the  system  variables  are  capacitor  voltages,  inductor  currents  or  other 
variables  determined  by  the  CAD  system.  Thus,  the  decision  to  use  a  CAD 
system  represents  a  trade-off  between  ease  of  system  formulation  and  some 
control  of  integration  error  by  judicious  choice  of  system  variables. 

The  state-variables  used  in  SCEPTRE  are  capacitor  voltage  and  inductor 
current.  Both  nonlinear  capacitors  and  nonlinear  inductors  are  present  in 
the  integrated  system.  The  nonlinear  inductors  arise  naturally  in  the 
generator  and  transformer  models,  the  nonlinear  capacitors  in  the  SCR  model. 
Further,  the  nonlinear  capacitors  in  the  SCR  model  are  strictly  monotonic 
so  that  computationally  the  choice  of  capacitor  voltage  as  the  state  variable 
is  desirable.  Also,  the  nonlinear  inductors  in  the  transformer  and  generator 
model  are  strictly  monotomic  so  that  inductor  flux  linkages  are  the  prefered 
state  variables. 

While  it  has  been  rigourously  proven  that  for  nonlinear  dynamic  networks 
containing  passive  elements  and  independent  sources  it  is  desirable  to  choose 
capacitor  charges  and  inductor  flux  linkages  as  state  variables,  no  one  has 
been  yet  able  to  prove  that  the  same  criterion  applies  for  networks  contain¬ 
ing  active  elements  and  controlled  sources,  although  it  would  seem  to  apply 
(Reference  11).  In  a  system  as  complex  as  the  one  in  this  project  it  is  not 
clear  what  the  trade-off  may  have  cost  in  terms  of  integration  error. 
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SECTION  VI 


NUMERICAL  METHODS  FOR  STIFF  DIFFERENTIAL  EQUATIONS 

6.1  INTRODUCTION 

In  this  section  suitable  numerical  methods  for  solving  stiff  differen¬ 
tial  equations  are  investigated.  Stiff  means  that  the  magnitude  of  the 
largest  eigenvalue  to  the  smallest  eigenvalue  of  the  system  differs  by  many 
orders  of  magnitude.  Most  standard  methods,  such  as  the  Runga  Kutta  method, 
or  the  forward  Euler  method,  are  not  well  suited  for  solving  stiff  differen¬ 
tial  equations  (Reference  15).  This  is  because  the  step  size  in  the  forward 
Euler  method,  for  example,  is  determined  by  the  smallest  time  constant,  but 
the  number  of  iterations  required  to  reach  the  steady  state  is  determined 
by  the  largest  time  constant.  Thus,  if  the  time  constants  are  widely  sepa¬ 
rated  a  small  step  size  would  necessarily  be  employed  and  a  great  deal  of 
computer  time  would  be  spent  on  determining  the  steady-state  solution 
(Reference  16). 

6.2  INTEGRATION  METHODS 

It  has  been  shown  (Reference  16)  that  the  widely  separated  time  constant 
problem  can  be  eased  by  using  implicit  integration  algorithms,  such  as  the 
backward  Euler  and  trapezoidal  algorithms.  Typically,  the  use  of  the 
trapezoidal  algorithm  is  avoided  because  of  the  "ringing  behavior"  it  ex¬ 
hibits  when  implemented  with  a  large  step  size.  This  is  undesirable  since 
it  may  erroneously  lead  the  designer  to  conclude  the  solution  to  the  system 
is  oscillatory  (Reference  16). 

Gear  (Reference  17)  has  proposed  a  family  of  methods  of  order  p+1  which 
have  been  shown  to  be  extremely  useful  for  solving  stiff  differential  equa¬ 
tions.  These  methods  are  based  on  the  implicit  formula 

yn+l  =  j0  clVl  *  dyrrtl 

For  p  =  0  the  backward  Euler  method  is  obtained.  The  Gear  methods  posses 
the  property  that  Re(hx)  <  -  a  for  some  small  positive  a;  that  is,  the 
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method  Is  stable  where  stability  Is  the  prime  consideration.  In  addition, 
when  Re{hx}  >  a  they  are  accurate  (Reference  17,  18).  In  general  there  Is 
a  tradeoff  between  accuracy  and  stability  for  any  numerical  Integration 
method.  Choosing  a  high-order  Gear  method,  for  example,  increases  the 
accuracy  (1e,  reduces  truncation  error)  at  the  expense  of  decreasing  the 
stability  (Reference  17).  The  maximum  order  must  be  limited  to  p+l*5 
because  of  stability  considerations  (Reference  17). 

Gear's  methods  can  be  implemented  in  a  fixed  or  variable  step  size 
mode,  or  in  a  variable  step  size  and  variable-order  mode.  While  the  latter 
is  difficult  to  implement  there  are  several  reasons  why  it  is  well  worth 
the  effort.  First  of  all,  since  the  computation  is  started  with  a  low- 
order  rule  and  small  step  size  the  method  is  self-starting.  The  order  and 
the  step  size  are  adaptive  in  the  sense  that  they  are  modified  at  each  time 
step  so  that  the  optimal  order  and  step  size  are  chosen,  subject  to  certain 
constraints.  This  leads  to  a  very  efficient  and  accurate  solution  over  the 
integration  interval  of  interest.  At  the  present  time,  stability  theory  in 
variable  step  size,  variable-order  methods  is  not  fully  comprehended  so 
that  decisions  for  determining  step  size  are  based  on  experimental  evidence 
(Reference  18). 

Another  advantage  of  using  Gear's  method  is  that  it  has  an  outstanding 
error  control  feature  (Reference  19).  Calahan  was  the  first  to  discuss  how 
to  implement  Gear's  method  in  a  nodal-based  analysis  program  (Reference  20). 

Gear's  method  is  a  vast  improvement  over  previous  techniques  because 
it  is  not  troubled  by  the  minimum  time  constant  problem.  However,  if  step 
size  is  changed  very  rapidly.  Gear's  method  could  become  numerically  un¬ 
stable.  To  overcome  this  deficiency  in  Gear's  method,  Brayton,  Gustafson 
and  Hachtel  (Reference  21)  have  proposed  yet  another  method  which  is  less 
likely  to  become  unstable  when  the  step-size  changes  rapidly,  as  demonstrated 
by  numerical  experiments.  In  the  case  where  the  step  size  is  fixed,  the 
method  of  Brayton  et  al.  can  be  shown  to  be  equivalent  to  Gear's  method 
Reference  21). 
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6.3  INTEGRATION  METHODS  OF  SCEPTRE 


SCEPTRE  offers  several  options  for  Integration  methods  including  Gear's 
(Reference  17)  method  in  variable  step  size.  SCEPTRE  should  be  appropriate 
for  this  project.  The  generator,  transformer  and  SCR  models  in  a  composite 
system  with  its  stiff  differential  equations  require  the  implicit  integra¬ 
tion  capability.  The  program  also  allows  specification  of  maximum  and  mini¬ 
mum  step  size  and  error  limits  can  be  specified  control  the  variable  step 
size  feature. 

However,  even  though  implicit  integration  allows  larger  step  size  than 
non-implicit  methods  there  is  no  doubt  that  a  step  size  appropriate  for  the 
short  time  constants  of  the  SCR  model  will  still  require  long  computer 
times  for  the  representation  of  the  periodic  swings  of  the  electromechanical 
equations  of  the  generator. 
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SECTION  VII 


SYSTEM  INTEGRATION 

The  system  integration  parts  of  this  project  did  not  have  any  specific 
effort  scheduled  for  the  period  of  this  report.  However,  in  making  the 
decision  to  use  SCEPTRE  as  a  modeling  and  simulation  program,  system  in¬ 
tegration  became  involved. 

The  original  plan  was  to  model  the  system  by  components  and  form  the 
system  by  integrating  the  components.  This  is  in  contrast  to  modeling  the 
system  by  elements  (R,  L,  C  and  sources)  and  integration  of  the  system  with 
network  equation  formulation. 

As  work  began  in  the  study  of  methods  for  solving  "stiff"  differential 
equations  several  problems  became  apparent.  First,  small  step  sizes  and 
long  computation  times  would  be  required  to  account  for  the  effects  of  both 
long  and  short  time  constants.  Secondly,  while  implicit  integration  methods 
for  "stiff"  differential  equations  allow  larger  step  sizes,  iteration  is 
required  for  convergence  to  a  solution  at  each  step.  Finally,  any  further 
iterations  that  may  be  required  to  find  the  equilibrium  point  for  the  com¬ 
ponent  interface  variables  make  the  component  approach  less  attractive. 

Three  alternatives  were  considered. 

1.  write  a  network  solving  program 

2.  use  CAD  program  SPICE  2 

3.  use  CAD  program  SCEPTRE 

The  first  alternative  has  the  advantages  of  a  much  simpler  and  smaller 
program  than  SPICE  2  or  SCEPTRE  bince  they  were  written  to  do  much  more  than 
a  simple  transient  network  analysis.  Also,  system  variables  could  be  used 
which  would  enhance  the  accuracy  of  the  nonlinear  systems.  (Note  Section  V) 
SPICE  2  and  SCEPTRE  normally  use  voltage  and  current  as  variables. 

The  disadvantage  is  an  unpredictable  amount  of  time  in  writing  and 
debugging  such  a  program. 

While  SPICE  2  is  well  recognized  as  a  good  computer  aided  design  program 
for  electronic  components,  it  has  limited  capability  in  representing  nonlinear 


time  varying  coefficients  such  as  the  inductance  of  the  generator  model. 
SCEPTRE  allows  definition  of  such  elements  through  functions  and  FORTRAN 
function  subprograms. 

SPICE  2  has  the  capability  of  defining  dependent  sources  as  polynomial 
functions  of  current  variables,  thus,  some  nonlinerities  can  be  represented. 
Possibly  the  time  varying  nonlinear  generator  inductances  could  be  represented 
in  SPICE  2  but  SCEPTRE  is  preferred  since  such  manipulation  is  not  required. 

Both  SPICE  2  and  SCEPTRE  have  implicit  integration  options  for  ’stiff" 
differential  equations.  Both  have  variable  step  size  with  some  external 
control  through  error  limit  specifications.  Also,  both  have  the  capability 
of  storing  electronic  component  models. 

In  SPICE  2  the  system  elements  are  supplied  to  the  program  and  the 
programs  integrates  the  elements  into  a  system  using  a  nodal  formulation. 

In  SCEPTRE  the  system  elements  are  supplied  to  the  program  in  much  the  same 
format  as  for  the  SPICE  2  program.  SCEPTRE  then  integrates  the  elements 
into  a  system  by  formulating  state  equations.  It  is  well  known  that  the 
nodal  formulation  is  generally  more  efficient. 

The  overriding  features  that  resulted  in  the  selection  of  SCEPTRE 
over  SPICE  2  was  the  limited  capability  of  SPICE  2  in  representing  nonlinear 
elements.  Both  generator  and  transformer  model  components  not  only  have 
nonlinear  elements  but  these  elements  are  subject  to  saturation  effects. 
SCEPTRE's  ability  to  define  elements  in  more  general  terms  and  to  use 
FORTRAN  function  subprograms  for  further  flexibility  were  necessary  in  the 
generator  and  transformer  model. 


SECTION  VIII 


GENERATOR  MODEL  RESULTS 


8.1  INTRODUCTION 

The  results  illustrated  in  this  section  represent  the  culmination  of 
many  experimental  computer  solutions.  At  each  step  in  the  development  of 
the  generator  model  two  types  of  runs  were  made.  First,  runs  were  made  to 
determine  appropriate  initial  conditions  so  that  subsequent  runs  would 
display  results  due  to  faults  and  not  due  to  mismatch  of  initial  conditions. 
These  initial  conditions  are  not  obvious  in  nonlinear  systems.  Secondly, 
runs  were  made  under  severe  fault  conditions  to  assure  that  the  model  would 
properly  represent  severe  transients  expected  in  later  simulations. 

The  steps  in  developing  the  model  included:  allow  the  speed  to  vary; 
add  torque  equations;  formulate  linear  model;  account  for  saturation  in 
terms  of  Equation  (13). 

In  accounting  for  saturation  for  the  first  step  the  and  similar 
terms  were  taken  as  zero  and  the  C  and  similar  terms  were  taken  as  unity. 

da 

This  gives  a  run  for  the  unsaturated  linear  model. 

The  next  run  was  made  where  the  Aa  and  similar  terms  were  taken  as 

a 

zero  but  the  C  and  similar  terms  were  computed  as  described  in  Appendix  d. 

ad 

This  gives  a  run  including  the  change  of  the  inductance  due  to  saturation. 
The  final  run  includes  the  variation  of  the  Aa  and  similar  terms.  This 

a 

will  include  the  effect  of  the  rate  of  change  of  inductance  with  current. 
These  terms  are  many  and  complex  and  a  successful  run  has  not  been  completed 

8.2  DISCUSSION  OF  RESULTS 

Figure  13  and  14  show  Cga  and  Cpp  curves  as  examples  of  the  saturation 
coefficients  used  in  the  example  of  the  subsequent  figures.  Data  for  the 
other  saturation  coefficients  are  included  in  Appendix  C. 

Figures  15  through  19  show  some  systems  variables  for  the  conditions: 


1.  one  cycle  of  steady  state  operation 


2.  a  symmetric  three  phase  short  circuit  is  applied 
at  the  stator  terminals  for  one  cycle 

3.  the  short  is  removed  and  system  load  is  restored. 

The  above  sequence  has  no  significance  except  that  it  represent 
severe  transient  under  load  and  the  model  responds  as  expected  in  all 
circuits. 

Figure  15  shows  phase  "a"  current. 

Figure  16  shows  the  main  field  current. 

Figure  17  shows  the  direct  axis  damper  circuit  current. 

Figure  18  shows  the  quadrature  axis  damper  circuit  current. 
Figure  19  shows  the  rotor  speed. 
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FIGURE  13.  SATURATION  MULTIPLIER,  C  VERSUS  Iv 
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2.  Implement  the  specific  AC  charging  configuration  to  be  used  in 


comparing  modeling  and  simulation  results  with  Air  Force  hardware  test 
results.  It  is  expected  that  this  will  involve  series  and  parallel  combina¬ 
tions  of  SCR's  to  accomodate  the  required  voltage  and  current  volves. 

Before  the  Phase  III  simulations  can  be  carried  out  the  specific  devices, 
the  number  in  series  and  the  number  in  parallel  and  the  exact  charging  circuit 
voltages,  currents  and  configuration  must  be  supplied  by  the  Air  Force. 

The  conversion  of  SPICE  2  models  to  SCEPTRE  models  can  be  construed 
to  be  part  of  Phase  II  task  3,  System  Integration. 
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FIGURE  16.  FIELD  CURRENT  Ip  FAULT  CONDITIONS  AS  IN  FIGURE  15. 
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FIGURE  17.  DIRECT  AXIS  DAMPER  CURRENT  FAULT  CONDITIONS  AS  IN  FIGURE  15 
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FIGURE  18.  QUADRATURE  AXIS  DAMPER  CURRENT  FAULT  CONDITIONS  AS  IN  FIGURE  15. 


APPENDIX  A 


THE  GENERATOR  MODEL 


The  purpose  of  this  appendix  is  to  establish  the  notation  and  address 
the  approximations  and  assumptions  in  the  generator  model. 

A. 1  EQUIVALENT  CIRCUIT  EQUATIONS 


Figure  20  shows  the  three  phase  generator  equivalent  circuit, 
vector  circuit  equations  are 


»  ’  RI  +  & 

where 


V  t-va*  vb.  vc,  Vp,  vD,  V/ 

*  ~  ^a*  ^b”  V  V*  ^D*  ^Q^ 

x  “  ^Xa*  xb*  Xc’  XF’  xD*  VT 


The 


(43) 

(44) 

(45) 

(46) 


(47) 


va,  vb,  vc  are  the  stator  phase  terminal  voltages 
Vp  is  the  rotor  field  terminal  voltage 
Vp  and  Vq  are  the  damper  winding  terminal  voltages 
ia,  ijj,  ic  are  the  stator  phase  currents 
ip  is  the  rotor  field  current 

ip  and  ip  are  the  direct  and  quadrature  damper  field  currents 

Xa’  xb*  xc’  XF*  XD’  XQ  are  the  c1rcuit  f1ux  1inka9es  with  the  subscrits 

referring  to  the  same  circuits  as  for  V  and  I. 
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In  terms  of  inductances  Equation  (43)  can  be  written 


V  =  RI  + 


3t 


(LI) 


with  the  flux  linkages  defined  by 


X 

a 

^aa 

Lab 

Lac 

LaF 

LaD 

LaQ 

>* 

cr 

Lba 

Lbb 

Lbc 

LbF 

LbD 

LbQ 

xc 

Lca 

Lcb 

Lcc 

LcF 

LcD 

LcQ 

XF 

LFa 

LFb 

LFc 

lff 

lfd 

lfq 

XD 

LDa 

LDb 

LDc 

ldf 

ldd 

ldq 

XQ 

_LQa 

LQb 

LQc 

lqf 

lqd 

lqq_ 
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(49) 


Some  approximations  and  assumptions  are  implied  in  the  choice  of  the 
inductance  terms  used.  These  involve  the  space  distribution  of  the  air  gap 
flux  and  the  resulting  variation  of  inductance  with  rotor  position.  The 
inductances  are  given  by  the  following: 


*"aa  Ls  +  Lm  Cos(2e) 

Lbt>  =  Ls  +  Sn  Cos  2(6  ‘  T> 

Lcc=  Ls  *  Lm  C«s  2<e  +T> 

Lab  '  -Ms  '  U  Cos  2<e  +  3-1  -  '■ba 

Lbc  *  -Hs  '  Lm  C°s  2le  -  T>  *  Lcb 

Lca  '  '"s  -  lm  c°s  2<e  +  X>  ’  Lca 
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=  L 


Fa 


LaF  =  Mp  Cos(e) 

LbF  =  Mp  Cos  (e  -  ^4  -  Lfb 

LcF  =  MF  Cos  <6  *  T1  ’  LFc 
La0  =  M„  Cos(e)  »  L 


Da 


LbD  s  M0  Cos  <e  '  T*  =  LDb 
LcD  •  Md  Cos  (e  +  -  LDc 

LaQ  *  Mq  Sin(e) 

L 


=  L 


Qa 


bQ  =  MQ  Sin  (0  “  3")  “  LQb 


LcQ 


Mq  Sin  (0  +  yO  =  L 


Qc 


LPP  =  Lf 

LDD  *  ld 

lqq  *  LQ 
lfd  =  mr  =  ldf 
lfq  =  0  =  lqf 
ldq  =  0  =  lqd 


(50) 


Figure  21  defines  the  rotor  angle,  e  used  in  the  inductance  expressions 

e(t)  =  ut  (51) 

A. 2  SATURATION  EFFECTS 

References  3,  5,  6  show  curves  of  the  inductance  variation  with  rotor 
position,  0,  similar  to  the  curves  of  Figure  22.  Kimbark  (Reference  5)  and 
Anderson  (Reference  6)  use  the  same  equation  form  as  used  here  in  Equa¬ 
tion  (50).  Smith  (Reference  3)  uses  an  L  _  of  the  form 


Laa  =  L*  +  Lm?  Cos(20)  +  ^  Cos(40) 


m2 


’m4 


(52) 
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FIGURE  21.  CONVENTION  FOR  DEFINING  ROTOR  POSITION 
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the  fourth  harmonic  being  due  to  a  third  harmonic  in  the  space  distrubution 
of  the  air  gap  flux.  Smith  further  states  that  this  fourth  harmonic  is  the 
highest  significant  harmonic  observed  over  many  tests. 

The  simpler  model  of  Equations  (50)  was  used  here  for  various  reasons. 

First,  the  inductance  model  of  Equations  (50)  is  a  generally  accepted  model 

and  it  is  simpler  than  Smith  (Reference  4).  This  becomes  significant  when 

the  non-linear  terms  of  tt  are  expanded. 

01 

Secondly,  it  is  clear  that  the  third  harmonic  in  the  flux  space 
distribution  is  evidence  of  the  peaks  of  the  wave  being  flattened.  This 
effect  is  due  to  saturation  effects  or  is  at  least  analogous  to  sturation 
effects  and  saturation  effects  in  our  model  are  being  accounted  for  by 
direct  measurements. 

In  this  work  saturation  is  taken  into  account  by  assuming  the  shape 
of  the  curve  of  inductance  variation  with  0  will  not  change  but  the  in¬ 
ductance  coefficients  in  the  equation  change.  For  example  the  saturated 

Laa  becomes  L  „ 
aa  aas 

Laas  “  Caa  Laa 

*  caa  Ls  *  caa  Lm  C°S<2|»  (53> 

Where  the  C  is  obtained  from  a  measured  curve  showing  the  variation 

aa 

of  L  as  a  function  of  a  net  excitation,  i„,  of  the  magnetic  circuit. 

ad  a 

It  is  clear  that  more  saturation  gives  a  larger  third  harmonic  in  the 
flux  space  distribution  and  hence  also  changes  the  shape  of  the  induc¬ 
tance  curves  of  Figure  22.  According  to  Smith  (Reference  4)  this  changes 
the  relative  amplitude  of  the  fourth  harmonic  compared  to  the  second. 

Smith  computes  his  saturated  L  as  follows. 

aas 


where 


'aas 


=  L 


si 


O+f) 


Caa  -  t-n  V  Caa  +  L„ 

- 9.  Cos(2e)  +  y  — — y - 9.  Cos(4e)  (54) 


Lj  is  the  direct  axis  inductance 
Lq  is  the  quadrature  axis  inductance 
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L 


90  180  270  360 

e  in  degrees 

Figure  22.  VARIATION  OF  INDUCTANCE  WITH  6 
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I 


K  is  the  decimal  fraction  of  third  harmonic  to  fundamental  flux 


Ld  Caa 


The  Equation  (54)  formula  assumes  saturation  effects  on  the  direct  axis 
and  no  saturative  effects  on  the  quadrature  axis  which  seems  reasonable  for 
a  salient  pole  machine.  Smith  does  not  refer  to  the  fact  that  K  is  a  func¬ 
tion  of  saturation. 

While  the  Smith  formulation  Equation  (54)  does  address  the  fact  that 
saturation  would  affect  the  shape  of  the  La.(e)  curve  it  is  quite  complex 
and  it  is  not  clear  that  it  would  give  better  results  that  Equation  (53) 
with  measured  C  a  since  it  does  include  several  assumptions  and  approxima- 

do 

tions. 


This  is  an  area  that  may  be  worthy  of  further  investigation  if  the  com¬ 
puter  simulation  results  do  not  match  observations  adequately. 


A. 3  THE  TORQUE  EQUATION 


The  equation  relating  electrical  circuit  quantities  to  mechanical  rota 
tion  quantities  is 


=  t  [Tm  -  Tp  -  B  If] 

2  J  m  e  dtJ 


(55) 


where 

J  is  the  polar  moment  of  inertia  of  the  rotating  parts 

B  represents  rotational  losses  and  may  include  viscous 
damping  and  possibly  core  losses  related  to  rotational 
speed. 


Tm  is  the  mechanical  torque  input  at  the  shaft 
Tg  is  the  so  called  electric  torque  developed. 


The  electrical  torque  developed  is  given  by  the  expression 

Te  =  [[I]T[|j][I]  (56) 

The  circuit  model  and  Equations  (55)  and  (56)  do  not  include  any  eddy 
current  or  hysteresis  loss  effects. 


These  losses  are  sometimes  considered  as  rotational  losses  and  the  B 
in  Equation  (55)  is  adjusted  to  include  them.  It  is  quite  common  to  neg¬ 
lect  both  the  core  losses  and  the  viscuous  damping  and  omit  the  B  term. 
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APPENDIX  B 


SATURATED  GENERATOR  MODEL 
B.l  SATURATED  GENERATOR  EQUATIONS 

Section  II  of  this  report  discusses  the  effects  of  saturation  on  the 
network  model.  Equations  (12)  and  (13)  illustrate  the  phase  "a"  equation 
development.  This  section  shows  more  detail  in  the  phase  "a"  equation 
and  also  lists  the  equations  for  the  other  phases. 

This  development  is  based  on  the  expansion  of  Equation  (43)  and  (48)  of 
Appendix  A.  From  these  we  write 

al '  A  (L1>  (57) 


d\  .  dl  .  ,dL«  T  (58) 

dt  '  L  dt  'dt'  1 

The  expansion  of  the  matrix  products  gives  as  the  first  row  of  Equation  (58) 


~ST  =  Laa  dT  +  Lab  W  +  L3C  dT  +  LaF  1ST 


din  dL  dL,h  dL 

Q  .  aa  4  A  ab  t  ,  ac  .. 

~st  +  ~sr  ja +  ~ar"  ^  +  ~at" 1 


'aQ  dt 


+  i  +^aif 


.  aF  . 

c +  ~sr  nF 


Saturation  effects  on  the  first  six  terms  on  the  right  hand  side  of 

Equation  (59)  are  as  shown  in  Appendix  A  Equation  (53). 

dLaa  ,  % 

The  development  of  the  on^derivative  term  is  shown  in  Equation  (12) 

The  complete  set  of  terms  for  is  now  listed 

d*a  T  »Lm1  3i  f  -  9Laal  9ib 

-ar  =  [NFa  cosle>  ’a  -H*[NFa  cos(e  ‘  X>  ’a  IT 


+  Nf#  cosfe  +  ^-J  iaTf 


+  Ncn  1 


FD  'a  TT~  TF 


-  [Nfa  *a  TT1  <*a  s1"^*  'b  si"<9  '  T>  *  1  c  sin(9  +  T» 

L  X  -a 


+  ,»  2Ca«Lm  S1"  29]lf 


|NFa  Cosl6)'b  TT^  ]  TT  +  |HFa  Cos(e  '  T1  'b  ~ f  ]  IT 

f  ....  2a,  .  3Lab  1  9lc  [  8Labl  s1F  .  L  ,  9Labl  9l0 
I  NFa  Cos(e  b  TT^  J  TT  |_'b  TT^J  IT  +  [NFD  *b  TT^J  IT 

["Fa  ’b  Ilf  «a  Si"l9)*  'b  Sin(®  -  T>  *  ’c  Sin(e  +  T» 

+  'b  2CaaSa  Si"<29  ’  T>]  II 

[»Fa  Cos  (e)  ic  ^1  ♦  [%a  costa  -  2f)  >c  ^ 

r  o  aL__i  av  r  3La_  3iP  r  3L  1  aio 

+  I  NFa  Cos(e  +  Q)  ic  -ft- J  -at-  +  [ic  -g^-j  +  |^FD  ic  iyj 


N  i  ^0 

TD  ’c  TT^  at 


-  |^NFa  ic  (ia  Sin(e)  +  ifa  Sin(e  -  +  ic  Sin(e  +  Q)) 

+  ic  2CaaLm  Si"<2e  *  T^]  H 
3L  c  3i  9  3L  F1  3ih 

+  Npa  Cos (0)  ip  -a^-  -atf  +  NFa  Cos(0  -  y)  nF  aix  J  IT 


["Fa  -.e  ♦  !?MF  &  ♦  [.F  3f]  £  *  [""  ^ 

[&L 

Npa  iF  -gy  (ia  sin(0)+  ib  Sin(e  -  Q)  +  ic  Sin(e  +  ^)) 

*  fF  CaFMaF  S1n(9)]  H 
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Cos(e)1D  -gT 


3ia  +  \h  rco<  *  2lt)  i  3LaDl  3lb 
IF  +  I  NFa(Cos  e_  T}  *0  TT^J  IT 


+  [“f, 

i  \n  Cn-fo  ,  2w\  .  3Laol  3ic  r  3LaD]  3iF  L  i  3LaDl  3lD 

+  I  HFa  Cos(e  +  T)  iD  -gY-J  -jf  +  |^1D  -gT^-  TF  pD  *D  TT^J  at 

*  [NFa  *D  (ia  Sin^+  ib  Sin(e  -  T}  +  ic  Sin(e  +  T)] 

Iae. 
at 


+  CaDMaD  Sin  0 


+  Npa  Cos(e)  ig  -g^. 


!!a 

at 


aL.nl  ai. 


NFa  Cos(0  -  £>  iQ  ^ 


+  |Npa(C0se+^)  iQy^ 


lie 

at 


aL 


'Q  ai 


aQ 


x  J 


!1f  +  L  i  3LaQ 
at  +  [nfd  7q  Ti 


x  J 


3iD 

TF 


[3L 

NFa  <q  irf  <(a  s,'"(9)t  <b  S1"<e  -  T>  +  'c  sin(0  +  T» 

'  'q  caQMaQ  Cos(8']  It 


dia  dib  di_  dic  din  di, 


'aa  dt  ab  dt 


Ul  U|p  Ulj.  U|Q 

+  Lac  ~F  +  LaF  St  +  LaD  ”F  +  LaQ  ST 


(60) 


To  make  the  expression  appear  less  formidable  the  following  terms  are  defined. 


aL. 

Aa  “  ^a  ~aT 


aL.K  aL.  aL  p  aL.n  <^L  n 

ja  ,  •  ab  ,  j  ac  ,  •  aF  ,  •  au  ,  •  ai| 

~  +  nb  TT"  +  ’c  ITT  +  F  ~dT~  nD  TT"  +  Q  “aTT  161 J 

XX  X  X  X  a 

Ba  =  Npa  [ia  Sin(e)+  ib  Sin(e  -  y)  +  ic  Sin(e  +  y)]  ,  (62) 

and  ca  =  2CaaLm  ^a  Sin(2e)+  *b  Sl'"(2e  "  T>  +  *c  Si"(e  +  T>]  (63) 
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The  voltage  generated  in  phase  "a"  can  now  be  written  more  compactly  as 
■ar  ■  [AaNFa  c°s<6>*  caa<Ls  +  L»,  c°8(28»]  af 

+  [flaNFa  Cos(e  '  T>  *  Caa(_Ms  +  Lm  Cos(2e  '  T1]  ar 

*  [AaNFa  Cos(9  +  X>  +  caa(""s  *  lm  Cos(2s  +  T1]  af 

*  [Aa  +  caFMaF  c°s<8)]  af  '  [AaNFD  *  CaD«aO  c°s<8>]  af 
'  CaQMaQ  S'n(o)af 

-  [BaAa  +  Ca  +  <VCaFHaF  +  *0  CaDMaD)  S1"<8>-  'q  CaQM.Q  c°8(8>]  j£  l64> 
For  the  other  generator  circuits  the  results  are  now  listed 

if  -  [AbNFa  Cos(e)-  CaaMs  +  CaaLm  Cos(20  '  T>]  if 

+  [AbNFa  Cos(e  -  Q)  +  CaaLs  +  CaaLm  Cos (2e  +  x}]  if 

+  [AbNFa  Cos<0  +  Q  '  CaaMs  ’  CaaLm  Cos(2e)]  if 

+  [Ab  +  CaFMaF  Cos(e  "  T-]  if 
+  [AbNFD  +  CaDMaD  Cos^0  "  T*]  if 
+  [CaQMaQ  Sin(e  "  X5]  if 

-  [AbB  +  2Cb  CaaLm  +  <CaFMaF  V  +  CaDMaD  V  Sin(0  ’  T 

"  CaQMaQ  ’Q  Cos(e  '  t]  3T  ^65) 


66 


where 


a  =  ,•  !hba  _|_  ,  ^bb  A  ,  3Lbc  A  ,  3LbF  L  ..  3LbD  .  .  3LbQ 
b  a  n,  b  31X  c  31X  ’F  31X  ’D-317  ’q-ST^ 

Cb  =  ia  Sin(2e  -  +  ib  Sin(2e  +  ^)  +  i£  Sin(2e) 


dx  r  -1  di 

lit  =  [AcNFa  Cos^'  caaMs  +  CaaLm  Cos(2e  +  12°)J  ST 

+  [AcNFa  Cos(0  *  120>  "  CaaMs  +  CaaLm  Cos^20)]  IT 

+  [AcNFa  c«(b  +  ’20)  -  0aaLs  +  caaLm  Cos<29  *  120)]  TT 

*  [Ac  +  CaFMaF  C°sle  +  120>]  7T 

+  [AcnFD  +  CaDMaD  Cos(e  +  12°)]  "ar 
+  [VaO  *  Si"(»  *  '20)]  IT? 


(66) 

(67) 


Where 


3L. 


3L, 


V  +  CaDMaDiD^Sin^e  +  12°) 

-  Maq  ’q  c»s(»  ♦  ’20)]  jjf 

(68) 

,Lcc(,  8LcF  +  j  SLcD  + 

C  TT7+  ’FT i~+  ’DT T"+  lQ 

>(«) 

X  X  X 

X 

,  Sin(2e}+  ic  Sin  (2e  -  120) 

(70) 

-SC  *  [AFNFa  Cos(e)*  CaFMaF  c°5<9)]  ~SC 

*  [AFAFa  Cos(8  '  ,2°)  +  caFMaF  Cos(e  '  12°)]  IT 


67 


[AFNFa  Cos(0  +  >20)  +  caFMaF  Cos(0  +  120)1 
[*  *  CFFLF]  £ 

[afnfd  +  cfdmr]  TF 

[BAF  -  CFCaFNaF]  H 


where 


Ac  =  i. 


aLFa  3LFb  3LFc  .  3LFF  3LFD 

+  'L  FT!  +  1  _  T3 +  lr  ~Z1 +  lr 


F  ■>  b  'c  ^7  F  D  nx 


dX, 


1xd  r  i  d\ 

dt  “  [aD  NFa  Cos  0  +  CaDMaD  Coste)J  -gt" 

h  [ad  NFa  Cos(0  -  120)  +  CaDW 

‘  [A0  NFa  Cos(e  +  12°)  +  CaDN 

■  [ad  +  cfdmr]  if 


VJ  1t 

- 

I  dih 

Cos(e  - 

120)  j 

|<nr 

1 

dic 

Cos(e  + 

120)  J 

C 

IT 

[nfdad  +  cddld]  "31 


ddldJ  IT 

iMaDCo]  ^ 


BAD  +  CaDnaOuD  I  3t 


3LDa  .  ,  3LDb 


3LDc  ,  ,  3LDF  .  ,  3LDD 


where  ^1^*1  '?  'uH 

CD  »  la  Sin  (e)  +  1b  S1n(e  -  120)  +  1c  Sin(e  +  120) 


(71) 

(72) 


CF  =  ia  Sin(e)+  ib  Sin(e  -  120)  +  ic  Sin(e  +  +  120)  (73) 


(74) 

(75) 

(76) 


68 


=  [VaQ  sin(0j]  at 


[caQ"aQ  Slnle 

-  120)1 

3ib 

at 

P 

-i 

3ie 

[v.(  si"(e 

+  120) 

. 

w 

at 

[Vaq]  I 

where 


C  =  i  Cos(e)  +  ih  Cos(e  -  120)  +  ic  Cos(e  +  120) 


1 


(77) 

(78) 
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APPENDIX  C 


GENERATOR  DATA  USED 

Considerable  difficulty  was  encountered  in  obtaining  data  for  machine 
model.  Data  were  needed  to  test  the  model  performance.  Some  data  were 
supplied  by  the  Air  Force  through  technical  publications  AFAPL-TR-75-87  and 
AFAPL-TR-77-31 .  In  both  publications  the  data  were  incomplete  even  for  a 
fairly  standard  model.  In  one  of  the  publications  the  data  purported  to 
have  been  used  was  inconsistant  in  the  sense  that  it  gave  mutual  inductance 
with  coefficients  of  coupling  that  were  greater  than  unity.  This  problem 
was  identified  after  spending  considerable  time  and  effort  and  having  several 
unsuccessful  computer  runs.  Further,  neither  of  the  above  reports  contained 
any  useful  saturation  data.  The  needed  data  is  not  available  currently. 

C.l  UNSATURATED  GENERATOR  DATA 

In  order  to  have  data  for  a  real  machine  and  move  forward  with  the  simu¬ 
lation,  data  were  taken  for  a  machine  from  Reference  6.  No  usable  satura¬ 
tion  data  were  available  from  this  reference.  Saturation  conditions  were 
simulated  using  trend  curves  from  Reference  3. 

The  machine  data  used  are  listed  below.  These  are  unsaturated  values. 

Balanced  three  phase  60  Hz  Generator 

Ls  =  4.152  mH 

L  =  0.074  mH 
m 

LF  =  2.189  H 
Ld  =  5.989  mH 
Lq  =  1.423  mH 
M$  =  2.076  mH 
MF  =  89.006  mH 
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L 


Md  =  4.721  mH 


Mq  =  2.269  mH 
Mr  =  0.108994  H 

Rated  Power  160  MVA 

Rated  Voltage  15  KV  (line  to  line) 

Inertia  Constant  H  =  2.37  sec. 

R,  =  0.001542  fl 

a 

RF  =  0.371  fi 

Rp  =  0.018421  fi 
Rq  =  0.018969  fl 

rated  field  excitation  voltage  =  375  volts 
C.2  SATURATION  DATA 

The  additional  data  required  for  the  model  are  saturation  data.  It  will 
be  necessary  to  measure  or  determine  in  some  other  manner  this  saturation 
data  for  the  specific  machine  being  modeled.  Measurement  methods  are  dis¬ 
cussed  briefly  in  subsequent  paragraphs. 

For  the  purpose  of  obtaining  numbers  for  the  machine  model  used  to  test 
the  computer  model,  curve  data^ere  used  from  Smith  and  Snider  (Reference  3,4). 
These  curves  were  normalized  and  applied  to  the  above  machine  data.  The 
premise  in  measurement  of  this  saturation  data  is  that  saturation  is  dependent 
on  some  net  equivalent  field  excitation 

*x  =  *F  +  NFD  +  NFa  j^a  Cos(0)  +  1b  Cos(9  "  X^ 

+  ic  Cos(e  +  y-)J 
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i 


\ 


(79) 


The  excitation  level  can  be  obtained  by  varying  only  ip  .  The  incremental 
inductances  are  obtained  by  perturbing  the  various  other  currents,  say  Aia  , 
for  incremental  L  .  Integrated  incremental  inductances  give  saturated 

aa 

inductances.  Thus  a  plot  of  Lfla  versus  ip  with  e  =  0  gives  appropriate 
saturation  data  for  the  "a"  phase  self  inductance.  Data  for  the  plot  of  Lfla 
vs  ip  is  normalized  and  put  into  an  algebraic  equation  form  using  a  least 
square  regression  analysis.  The  following  form  is  obtained  from  the  Laa 
versus  ip  data 

Caa  =  a10  +  all  ix  +  a12  +  a13  ix  t80) 

Saturation  data  must  be  determined  for  Lfla  ,  Lflp  ,  ,  Lap  ,  Lpp  .  From 

these  the  regression  coefficients  for  the  C's  are  found:  caa »  Cpp »  ldd * 

CaF  ’  CAD  *  CaQ  *  CFD  * 

Table  2  lists  the  coefficients  used  in  the  simulation.  These  were 
obtained  by  reading  data  from  curves  given  by  Smith  and  Snider  (Reference  3) 
and  are  not  necessarily  correct  or  typical  for  the  machine  used  here.  The 
ballistic  method  described  in  Snider  and  Smith  (Reference  3,4)  permits  the 
measurement  of  saturated  incremental  inductances.  These  are  integrated  into 
the  saturated  inductance  curves. 

C.3  MEASUREMENT  OF  THE  UNSATURATED  INDUCTANCES 

1.  L  and  L 
s  m 

With  the  rotor  locked  into  position  (e  =  0°),  measure  the  induc¬ 
tance  of  phase  A  of  the  stator.  The  measured  inductance  is  L$  +  Lm  .  Rotate 
the  rotor  and  lock  into  position  (e  =  45°).  Measure  the  inductance  of  phase  A 
of  the  stator.  The  measured  inductance  is  L$  . 

Rotate  the  rotor  and  lock  into  position  (e  =  90°).  Measure  the 
inductance  of  phase  A  of  the  stator.  The  measured  value  is  L$  -  Lm  . 

Compare  all  measured  values. 
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TABLE  2.  SATURATION  COEFFICIENTS 


ajo 

ajl 

aj2 

aa 

0.98367 

-2.9246 x 10"5 

-1 .0771  x 10"7 

FF 

1.0088 

-2.9867x  1  O'5 

-1.0519  xlO"7 

OD 

0.97827 

-2.1569  xlO'5 

-7.9586x  10"® 

AF 

1.0000 

-6.9202 x 10'5 

-8.4819x10"® 

AD 

0.9886 

-8.0120  XlO'5 

-6.1443x10"® 

AQ 

0.99419 

-2.6322  xlO"5 

-5.5511  xlO"® 

DF 

1.0007 

-6.9355  xlO"5 

-6.6663  x  10"® 

'xx 


=  ajO  +  ajl 


+  aj2 


.2  ,  .3 

x  j3  x 


aj3 

2.0394  xlO"11 
1.8699 x  10"11 
1 .5930  x  10"11 
1.6303  xlO"11 
1.2983  xlO'11 
1 .0182  x  10-1 1 
1.3054  xlO-11 
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2.  I_s  and  M$ 

With  the  rotor  locked  into  position  (e  =  105°),  connect  a  sinu¬ 
soidal  voltage  source  between  phase  A  and  neutral  of  the  stator.  Measure 
the  RMS  current  flow  in  phase  A.  Measure  the  resultant  RMS  voltage  induced 
between  phase  B  and  neutral.  The  mutual  inductance  between  phase  A  and 
Phase  B  is  computed  as 


„  WRMS> 

"Ms  '  »l"(RMS)  ’ 

a 


where  u>  is  the  radian  frequency  of  the  phase  A  current.  The  quantity  M$ 
whould  be  positive  and  should  be  smaller  than  L$  . 

With  the  rotor  locked  into  position  (e  =  60°),  measure  and  I  as 
before.  The  result  yields  Lm  _  Ms  •  With  the  rotor  locked  into  position 
(e  =  150°),  measure  Vb  N  and  Ia  as  before.  Calculation  results  in  -M$  -  Lm  . 
Compare  values  of  Mg  . 


Measure  the  field  coil  inductance  directly. 


voltage 
current 
phase  A 


With  the  rotor  locked  into  position  (e  =  0°),  connect  a  sinusoidal 
source  between  the  field  terminal  connections.  Measure  the  RMS 
flow  in  the  field  coil.  Measure  the  RMS  voltage  induced  between 
and  neutral  of  the  stator.  Calculate  MaF- 


Va.N(RMS) 
u)Ip(RMS)  ’ 


5. 


R 


0 


As  measured  from  the  field  coil  terminals. 


Vp(RMS) 
ZF  =  Ip (RMS) 


R  +  JujL 
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“2MqF 

Zp  =  Rf  +  JwLp  +  R— +  juld 

-  I R  U>2RDMDF  \  “3MDFLD 
"  '  RF  +  n2  I  2,2  |  "  1  „2  .  72,  2 


R5  +  -  ld 


RB  +  “  L5> 


=  Rp  +  Rg  +  J(o(Lp  +  Lg) 

Measure  the  field  coil  resistance  directly. 

Ra  =  R  -  Re 
e  r 

Le  ‘  L  •  LF 

Assume  a  value  of  damper  resistance  Rg 

Re(R§  *  “2<-o> 

DF  J  ~2^ 


m; 


6. 


LeRD 


Note  that  MR  and  MpD  are  the  same  parameter. 
MaQ»  H 

As  measured  from  phase  A  terminals, 

V  (RMS) 


Za  =  I  (RMS) 

a 


2m2 

"  V 


=  R  +  JtoL  ,  Z_  =  R.  +  Ju>L  +  D  .  , 
a  a  a  Rq  +  JuLq 


All  measurements  are  taken  with  the  rotor  locked  into  position  (e 


Measure  the  phase  A  resistance  directly. 
Rfi  =  R  -  Ra 


Le  -  L  -  La 


=  90°) 
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APPENDIX  D 


THE  SCEPTRE  GENERATOR  PROGRAM 


D.l  INTRODUCTION 

The  following  pages  contain  the  listing  of  the  SCEPTRE  program  for  the 
generator  model.  The  program  is  an  implementation  of  the  equations  of 
Appendix  B.  The  circuit  diagram  of  phase  "a"  is  shown  in  Figure  1. 

It  should  be  emphasized  that  this  is  a  working  program  and  not  a 
finished  product.  Somewhat  more  efficient  use  of  defined  parameters  will 
be  necessary  in  a  final  program.  Also,  the  particular  outputs  used  are 
for  program  analysis  purposes  and  not  necessarily  important  in  system 
simulation. 

The  first  four  pages  contain  the  circuit  model  and  the  run  controls. 
The  last  two  pages  contain  FORTRAN  subroutines  used  to  compute  the  torque 
relationships  and  the  nonlinear  saturation  effects. 


I 

a 


M 

;  i 


t'i 


0.2  Circuit  Model 


CIRCUIT  DESCRIPTION 
ELEMENTS 

m*i  3  a-  := pe  v 

PAL#  1  31-  j  =  T1 
RRL# 1 32-C-T1 
RCL# 1 33-0=T1 
RA#1  - 1 1  =  PRP 
Pb#2-12=PRP 
RC#3-13=PRP 

P33#  1  3-23=  EQUATION  3  (-  1  .  #PLM#PM  2  P3  #  PC  A.  A) 

R22#  1  2-22=  EQUATION  3  (-  1  .  #PLM  #P  2  P  3/ PC  AA  > 

R11  #  1 1-21=  EQUATION  3  C- 1  .  #PL M /0 . 0D3 ,PC A  A ) 

RF#4-14=  PRF 
RD,5-15=  PRD 
RQ#6-16=PRQ 

E12#21-31=EQUATION  3  (IL 22 , PL K#PM 2 P3 # PC A  A ) 
E13#31-41=EQUAT10N  3  ( IL33#PLN#P2 P3# PCAA ) 
E14#41-51=EQUATION  4  (IL44#PMAF#Q.0D(j#PCAF) 

E  1 5  #-  51-61=EQUATION  4  ( IL  55,  PM  A  D#f*.  OD  C#PC  AD) 

£16#  61-71=EQl)ATION  5  (  IL  66# PM  AG#  fl .  OD  #  -1  .  DO#  P  C AQ) 
E18#  71-91=EQU  ATION  11  (P  C)  F A  A  #  FN  F  A  #0L22#  F  M2 P3 > 

E19#  91-101  =EQUATION  11  (  P C  OF  A  A  ,PN  f  A  ,  D L3 3  , P2 P 3  > 

E101 #101 -1 11 =EQU AT ION  12  ( -1 . GDO #P CD F AA # 1 .OD f ,DL 4 4 ) 
E 1 1 1 #111-121 =EQUATION  12  (-1 . ( DO, PCOF AA# PNFD ,DL5 5) 
F121  #121-131  ^EQUATION  12  (  1  .  -JD^PCCF  AA#  PC  OF  B  A,PS  ) 
E21#22-32=EOUATION  3  ( IL  1 1 #PL M#PM 2P3 # PC A  A > 

F23#32-42  =  EWUATI0N  3  ( IL 33 #PLM #0 . ODD #PC A A  ) 

E24#42-52  =  EQUATION  4  (  IL44#PMAF#P12P3#PCAF) 

E 25 #52-62 =EQUATION  4  ( IL 55 # PM  AD# PM 2P3# PC  A D> 
E26#62~72=FQUATION  5  (  IL  6  6#PM  AO#  PM2P  3,-1  .  f>0#  PC  AQ  ) 
E28#72-92=EQUATI0N  11  (P O F A F #PN F A # DL 1 1 # 0 . ODD)  ‘ 
E29#92-1  02=EQUAT1  ON  11  (  PC  OF  AF:#PnF  A,DL33#P2P3) 
E102#1C2-1 12=E0UATI0N  12  ( -1  .  000  #  PCO  F  AP#  1  .CJD  C#  DL4  4) 
E112#112-122=ECUATI0N  12  (-1  .  f  DO  #  PC  3  F  AB  #  f’L  FD  #DL  5  5) 
F122  #122-132=EQUAT ION  12  (1.0D'#PCCF AF#PCOFBA#PS) 
E31#23-33=EQUATION  3  (IL  1 1 #PL M #P2 "3# PCAA ) 

E32#33-43  =  EQUATION  3  (  IL  22  #P  L  A.  #0  .  I  DO  #PC  A  A  ) 

E34# 43-53=EQUATION  4  (IL 44 #PM. AF#P 2 F 3 #PC A F > 

E35# 53-63=EQUATION  4  (IL 55#PM,AD#P2P3#FCAD> 
E36#63-77=EQUA TION  5  ( IL 66# PM AQ# P 2 P3 # -1 . DC# P C AQ ) 
F38,73-93=EQUATION  11  (P CO F A C # FN M ,DL1 1 # 0. DDO) 
L39#93-1 03=EGUATI ON  11  ( PCOF A C #PN F A , 0L22 # PM2 P3 > 
t103#  103-1  13=LQUATI0N  12  ( -1  .  ODO  #  r  CD  F  AC  #  1  .  OD {  #  DL  44  ) 
E113#113-123=C«UATI0N  12  ( -1  .  ODC-  #  =  C  0  F  AC  #  P  »F  D  #DL5  5) 
5123#  123-133  =  FQUATION  12  (  1  .  0  D  ’  #P  C  OF  A  C#P  CLFB  A#PS  ) 
E41#14-24=EQUATION  4  (  I L  1 1  # P  K  A  F# T  .  OD  #PC  A  F  ) 

E42# 24 -3 4= EQUATION  4  ( IL 2 2 #PN A F#P 1 2P 3 #PC AF ) 

C43#  34-44  =  EQUATION  4  (  IL  33#  PM  A  F#  P  2  P3  /  PC  A  F  ) 

E4S# 44-94  =  EQU ATION  11  (P C OF  A F #FWF A ,D L 1 1 # C . OD 0) 
C49#94-1l4=£UUATI0N  11  ( PC ~F A F #PN F A # D12Z # PM2P3 > 

E 104  # 1 04 -1 14  =  FOUA  T  ION  11  ( PC  0 F AF # PNF A #DL 3 3# P ?P3 ) 
E114#114-124=EQU  ATION  12  ( -1  .  ODC  #  PC  0  F  AF  #  PNF  D  #DL  5  5) 
Fl24, 124-1  34  =  FGUATION  12  ( 1  .  v'D’1#  PC  OF  A  F#  PC  OF  B  A#  F'S  ) 
651,1 5-25=EQUATI0N  4  ( IL 1 1 #PM, A D# n . #p C AD ) 
t52#  2c-3  5  =  EQUATION  4  ( IL  22  #PM  A  C#  nM  2P  3  #PC  A 10 
E  53  #  3  5  -4  r  =E  Q  U  A  T I  0  N  4  (IL33#FNAP#P2P3#PCAD) 
E58#45-95  =  EQU*TI0N  11  (P  COFAD  #FNF  »,DL11  #C*.0DO> 

E59# 95-1  ■  5=EQUATI0N  11  ( PCOF * D #PN F A, DL 22 , F M 2 P3 ) 
C105# 1C5-1 15  =  FGU AT ION  11  ( PC  0 F AD # PNF A#DL 3 3#P ?P 3) 
f115#115-125=EQU ATION  12  (-1 . 0  CO# PCO F AD# 1 . ODC#DL 4  4) 
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D.2  Circuit  Model  (Continued) 

E125,125-C=EQUAT1PN  12  (  1  .  ?D  .  ,PC0  f  AD  #  PCO  FF,  A,PS  ) 

E61# 1 6-26=EQUA TION  5  t  IL1 1  *P« AQ,0 .D, -1 .ODf,P CAQ > 

E62, 26-36=EQUA TION  5  <T(  "'2,',,,AO*PM2P3,-1.I}D0»PCAQ) 
E63#  36-0  =EQUATI0N  5  ( 1L 3 3, PM  A C, P 2P3 a -1 . 000, PC AQ ) 

RAD, 1 31-0  =  1 . DD4 

ROD,  1  32-0  =  1  .004 

RCD,  133-0  =  1  .004 

M12»L11-L22=PM12 

M13, L11-L33=PM13 

M14,L11-L44=PM14 

M15#L11-L55=PM15 

M16*L11-L66=PM6 

K23,LZ2-L33=PM23 

M24,  L22-L44=PM24 

M25,L27-155=PM25 

M26,L  22-L66=PM26 

M34, L  33-L4  4=PM34 

M35,L33-L55-PM35 

M36,L33-L66=PM36 

M45, L44-L55=PM45 

L11, 0-1  =  PL1 1 

L  22, 0-2  =PL22 

L33, C -3  - PL3  3 

L44, C -4  - PL4 4 

L55,C-5  =  PL5  5 

L66, 0-6  =PL66 

DEFINED  PARAMETERS 

P2P 3  =  2. 0943951 024  DC 

PM2P3=-2 . 094 395102 4D0 

PRQ= 1 8  .969D-3 

PRD=1f.421D-3 

PRP= 1 . 5420-3 

rRF  =  0.61  DO 

PLS=4 .470050-3 

PLM  =0.07  43  3D -3 

PLFF  =2.1 89D0 

PLDD  =  5 .93  9D-3 

PLQ3  =  1.4230-3 

PMS=1 . 75945D-3 

PMR  =  0.09  901383  DO 

PKAF  =89. 0060-3 

PMAQ=1 .96904D-3 

PMA5=4.22Q98D-3 

PRL=1 .75  78100 

P FV=  3  75.  ODO 

PCN  A  A  =  0. 9836  7  DO 

PA11 =-2.92460-5 

PA12=-1.0771D-7 

PA1  3  =  2  .  03  94  D  - 1 1 

PA1  4  =  0.  DO 

PCNAf  =1. OOOODO 

PA21  =-6.92120-5 

PA22=-8.4719D-8 

PA2  3  =  1. 63030-11 

PA24=  0.00 

PCN  A  0  =  0.  9886D0 

PA31  =-P. 0120D-5 

PA32= -6. 1 443C-6 

P A33  =  1.29330-11 

PA3  4  =  0.00 
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D.2  Circuit  Model  (Continued) 


PCNF o=1. DO 
PA41 =-6.93550-5 
PA42=-6.66630-F 
PA43  =  1 .3054D-1 1 
PA44  =  0.00 
PCN 00=0.9782700 
PAS  1  =-2.15690-5 
PA52  =-?.  95d6D-8 
PA53  =  1.59300-11 
P  A54  =  0.00 
PCNF F  =1. GO'S  DO 
PA61  = -2 . 93  67  D -5 
PA62=-1 .05190-7 
PA 63  =  1 . 8 69 90  - 1  1 
P  A64  =  0 . 0  P- 
PCN  AQ  =n. 9941 900 
PA71 =-2.fc322D-5 
PA72  =  -5.5S11  D-8 
PA73  =  1 .016  20-11 
P  A  7  4  =  0 . 0  0 
rfJFA=4.06  6D-4 
F’NFO  =  4 .9  30 07  L-2 

PCEa  =  FIEti(IL11,IL2  2,  1L33,IL44,IL55,PNFA,FNFD,PT) 

PCA A  =  FPOLY( PA11, PA12/PA1 3 /PA  14/PCta  ,3.03/PCNAA) 

PCF  F  =  FP0LY(PA61,PA6  2,rAt.3,PA64,PCEa,3.D3,PCNFF) 

PCDD  =  FPCLY<PA51,FA'2,FA53,PASA,PCEa,3.D3,PCN'DD) 

FCAF=  FP0LY<PA21,FA2?,PA23,PA  24/ P  C  f  Q  ,  3.  D  3  ,P  C  9  A  F  > 

PCAD  =  FP0LY(PA31/PA32/PA33/FA34/I CE0/3.03/PCNAD) 

PCA 3=  FP0LY<PA71/PA72/PA73/PA74/PCEa^3.D3/PCNF0) 

PC  F  3  =  FPOLY(PA41,PA4  2/PA^3#PA4  4/PCfO*3.03/PCMFO) 

F  OA  A  =  FPOLYl  (  P  Al  1  ,  FAl  2/P  M’A/P  A  H  ,  »  C  t  0 /3 . 03> 

P  OF  F  =  FPOLYl  (FA61/FA62/P  A63,PA'4,PCEQ,3.D3) 

PDDD=  FPOLYl ( PA51 / FA 52/P A53, PA54/PCE 0/3. 03) 

POA  F  =  FPOLYl  <  PA21  /  l;A??,P  A2?/  PA  24/  °CE  0/3.03) 

PDAO  =  FPOLYl  (FA31,  PA  3  2/ P  A  3  3,  F  A  34,  P  C  E  0,3.0  3) 

PDA3  =  FPPLY1(PA71,rA72/PA73,PA74/PCEQ/3.03) 

PDFD  =  FPOLYl  (PA41, F A 4 ? ,P A4 3, P  /  44 , PCE 0 ,3 . D 3) 

PCOF  R A  =  FS ItP  (PNF  A  ,IL  1  1/  IL2?  /  I L3  3/PT  ) 

PCOF  A  A=  FfAESSI  CIL1  1/1L22/IL33/IL44/I  L  55 /  I  L66  ,PD  A  A  ,  P  0  A  F  « PD  AD,  P  DAQ 

F’LS,PLK,F1S,PT,?.DDC,PP'AF,F''!AD,P«IAa,D'  ‘ 
PCOFAB=  FWES S 1 ( IL 1 1/ U 22/IL33/IL44/I L 55/ IL66/PDAA,P0AF,P0AD/P0A0 

FLS/PLK/P«IS/rT/PM2P3/PMAF,P'A  AD/P*1  S8 , 2)' - 

PCOF  A  C=  FMES  SI  ( IL1 1/ II 22 /IL3  3/IL4  4/I L55/ 1L66/P0A A/POAF /PD AO/F  OAB 

PLS/PLM,r'1S/PT/P2P3/P'1  AF/PMf  0/P«IAQ/3)~ — ’ 

PCOF  A  F  =  F*!FSS2(IL11,I122,IL33,!L44,IL55,  PDAF/PDAF/PDFD, 

PT/  PF’A  F,  PL  F  F,P*AR) 

PCOF  AD=  FMES  S2CIL  1  1/1  L22/IL3  3/IL4  4/I  L55/P0A.0,PDF  D/PDDP, 

F'T/PiXA  D/PMR  ,PL00  ) 

FM12  =  EQUATI0N  2  ( F MS , F LM , P C A  A  ,  PI 2 P 3 ) 

P*13=ECUATI0N  2  <f'9S/PLK,°CAA,P2P3) 

PP14=EOUATION  9  (1  .ODD, PKAF/PCAF, 0.000) 

Ph15=EQUATI0N  V  (1  .CDC/PMAO, PCAD, 0.30(3) 

PM1  5  =  E  Qll  AT  ION  8  (  1  .OD C ,PM AC, G  . G/PC A3 ) 

PP23=ECUATION  2  (PMS/FLM/PCAA/O.OOO) 

Py24  =  E  QUA  T  ION  9  <  1  . 00  2,  P F  A  F  ,  P  C4F  ,  P *2  P  3) 

P"25  =  EOUATION  9  (  1  .00  f  ,PP  AO  ,  P  C  AO  ,  Pf.2  P  3) 

PM25 = EQUATION  8  (  1  .0  DO,P  .1  A  0,  PK  2P  3,  PC  AQ) 

PP34  =  EQUATION  9  (1 .QDC/PK4F/PCAF,P?P3) 

PM35  =  EQUATI0N  9  ( 1  .GO C ,PM AD / P C AD / P2 P 3  ) 

P"35=EQUATION  8  (  1 .GO G,P K AQ, P 2 P3, PC  A  Q  ) 


0.2  Circuit  Model  (Continued) 


PV45=X13(PNR*FCFD) 

PU1  =FQUATION  1  (PLS/PLM/T.COC/PCAA,  PCOF  AA/PNFA) 

PL22  =  EQL'  *T  ION  1  (  P  LS  /  PLM /P  2P  3 /PC  A  A /P  C  OF  A  B /PN  F  A  ) 

PL33=EQUATION  1  <  PLS/ TLM /f>M2P  3/PC  AA/ PCOF  AC /PNF  A) 

PL44  =  X11(PLFF*PCF  F+PCUFAF) 

PL55=  Xl  2  (  PL  D  D*PC  D  0  ♦  PC  OF  AO  * FNF D > 

PL65  =  PLQQ 

DPS=X1<F02THT(IL11/II?2/I133/VRAD/VRPD/VRCD/PRP/PS>> 

DPT  =  X2( PS) 

PS  =3  76. 9980  0 
PT=  2 . 489490  0 
FUNCTIONS 

EQUATION  1  (A/P/C  / D/E/F)=(D*(A+B*DCDS(2.  Ci  DQ*PT  +  C))+  F * E* 0  COS ( PT-C) ) 
EQUATION  2  (A/B/C/D)=(C* (-A  ♦  P*D C OS < 2 . 00 C* P T+ D ) ) ) 

EQUATION  3  (A/B/C/D)  =  (?.D3*A*B*D*PS*DSIN<2.DC*PT*C)  > 

EQUATION  4  C  A / B/C / C)  =  (A*u*D*PS*DSIN( PT+C )  ) 

EQUATION  5  (A/B/C/D/f)=(£*D*A*D*PS*DCOS(PT+C>) 

EQUATION  7  (A/P/C) =  (  2 ,*A*B*DSIN(2 ,*PT  PC) ) 

EQUATION  8  (A/B/C/P)  =(A*6*D*0SIN(PT*C)) 

EQUATION  9  ( A / P/ C / 0) = < A* 0 *C * D COS ( P T* D  )  ) 

EQUATION  11  (A/B/C/D)  =  (-1.'>n''-*A*B*C*  DCOS ( PT+D) ) 

EQUATION  12  (A /P/C /D  )=  (A  *B*C  *D) 

TABLE  1 

O.ODO/1. 75  781  0  0/ j. Cl  6 666 0  0/1  .  757 81  DO / C.D 1 666 7DO / 0. 5  DG/ 0. 03333  DO/ C .50 0 

n  .03  3  340  n/1 . 7  5  781  D  b/  i'  .  130  0/1  .7  578100 

OUTPUTS 

I  Li  1  /  IL44/IL  55/IL66/PS 
PLOT  INTERVAL  =1.0-8 
INITIAL  CONDITIONS 
IL11 =6967.4400 
IL22=-3483.72D0 
I L33  =  -34  E  3 . 72  DO 
IL44  =  616 .31D0 
IL55=C.0DC 
I L65  =  O.CDO 
RUN  CONTROLS 

INTEGRATION  R OUT  I N E= I  PPL  1  C I T 
MINIMUM  AESOLUTE  ERROR  =  1.0-4 
MINIMUM  STEP  SIZE  =1.0-24 
S  TO  3  T  IN  E  =  0. 05  DO 
END 
/* 
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D.3  Fortran  Subroutines 


//FORT. FOR TSRC  DO  * 

DOUBLE  PRECISION'  FUNCTION  F  D  ?T  H  T  <  I  L  1 1 »  IL  2  2 ,1  L  33  ,  VR  A  D  ,  V  R  UD» 

1  VRCD,PPP,PS) 

IMPLICIT  RE AL*8 < A-h*0-Z > 

REAL  *8  IL1 1  *IL22,I  L33 

C  THIS  SUBPROGRAM  COMPUTES  ROTOR  ANGLE  ACCELERATION.  THE  MECHANICAL  AND 
C  ELECTRICAL  POWERS  ARE  CALCULATED  IN  P.U. 

DATA  SB3,PM,HCON,OMEGAR/160.OD6,0.8D0,2.37DD,376.998D0/ 

PE  1  =  IL 1 T  *VR  A  0  ♦  IL  22*tf  RBD  ♦  IL33*VRCD 

PE  =  PEI  ♦  PRP*UL11*IL11  ♦  IL  22*1 L22  *  IL33*IL33> 

PE  PU  =  PE/SB  3 

OMEGAU  =  PS/OMEGAR 

DENOM  =  2.0DUH  CON  *3M£  G  AU 

DIFF=PK-PEPU 

FD2THT  =  DIF  F/DENOM 

RE  TURN 

END 

DOUBLE  PRECISION  FUNCTION  F  I  £Q  (  ILl  1,1  L  22,  I  L  33 ,1  L44,  IL5  5  , 

1  PNF A,PNFD,PT> 

IMPLICIT  PE  AL*8 (A-H»I*0 -Z) 

C  THIS  SUBR 0 G pAM  COMPUTES  THE  EQUIVALENT  CURRENT  ALONG  THE  ROTOR  AXIS. 
C  REFER  TO  PAPFR  BV  SMITH  AND  SNIDER  FOR  THEORECTICAL  FORMULATION. 

A  =  DCOS(PT) 

B  =  DC0S(PT_2. 0943951  02393  DO) 

C  =  DCOS (PT  *2.094395102  393 DO) 

ST  £  PI  =  < I L 1 1  *  A  +  I L  22  *  B  ♦  IL33*C)*2.D0/3.D0 
STEP2  =  PNF  D  *1 L  5  5  ♦  PNFA*STEP1 
FIEQ=  IL44  ♦  STEP2 
FIFOs  DAhS(FIEO) 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  F  POL  Y  ( A 1  /■  A2  /•  A3  ,  A  4  r  C  EO  r  PL  I  M  C  0  NS  ) 
IMPLICIT  REAL*R(A-Z) 

IF  (CEO  .GT.  FL 1 M ) C  EQ  =PL  1M 

CEQ2=CEQ*CEQ 

CEQ3=C  EQ  *CE  Q  2 

CE04  =C  EQ  *CE  Q  3 

FPOLY  =  CONS  ♦  A  1 *C  E  Q  ♦  A  2 *C E 3 2  ♦  A3*CEQ3  *  A4*CCQ4 

RETURN 

END 

DOUBLE  PRECISION  FUNCTION  F  P  OL  Y  1  (  A  1  ,  A  2  ,  A  3,  A  4,  C  E  Q,  PL  I M) 
IMPLICIT  REAL*? (A-Z) 

CE02=CEQ*CEQi 

CEQ3=CEQ*CEQ2 

CEQ4=CEQ*CEQ3 

F  P  OL  Yl =  A1  ♦  2.P0*A2*CEC  ♦  3.Dr*A3*CEQ?  ♦  4.D0*A4*CEO3 
IF  ( C EQ  .GT.  FL IM)F  POLY  1 =n.D0 
RE  TURN 
END 

/* 

// 
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D.3  Fortran  Subroutines  (Continued) 


DOUBLE  PRECISION  FUNCTION  F  *1 E  S  S  1  (  I  L  1 1,1  L  22  ,  I  L  33,  I  L  4  4  ,  1  L  55,  I  L  66  , 

1  XLAA,XLAf  ,XL4D,XL4e 

2  PLS,PLM,PMS,THT,PHASE,PMAF,PMAD,PMAQ,NFLAG> 

IMPLICIT  REAL*P(A-H,I,0-Z) 

C  THIS  SUBPROGRAM  COMPUTES  COEFFICIENTS  OF  VOLTAGES  SOURCES  PUT  IN  TO 
C  TAKE  INTO  ACCOUNT  THE  EFFECTS  OF  SATURATION. 

GO  TO  (  10,20,30),NFLAG 
10  XI  =PLS 
X2  =  -PM  S 
X3=-PMS 
GO  TO  40 
20  X 1  =  - pM S 
X2=  PL  S 

X  3  =  -PM  S 
GO  TO  4  C 

3  C  X1=-PMS 
>2  =-PMS 
X3=PLS 

40  DL  A  A  =  XLAA*(Xl  +  PL  A  *D  CO  S  (  2  .  OD  C  *TH  T  +PH  AS  E  )  ) 

DL  Bb  =  X  L  AA*  (  X2  ♦  P  L  A*  0  C  C  S  C  2  .  OD  0*T  H  T -2 . 0^4  3  951  0 2  393  D  C  *P  H  A  SE  )  ) 

DLCC=  XLA.A*<*3  +  PLM*DCOS<2  .GD  C:*THT  +  2. 09439 51 02S93D0+ PHASE >> 

OL  A  F  =  XLAF*P!',AF*DCOS(THT*PHASE) 

DL  AD  =  XLAD*PNAD*DCOS(THT+PHASE > 

DL  AQ  =  XLAQ*PMAQ*DS IN  (THT+PHASE  ) 

FMESS1=  IL 1  1  *D  L  A  A  +  IL2?*DL3B  +  IL33*DLCC  ♦  IL44*DLAF  ♦  I L  5  5*D  LAD 
&  +  IL66*DLAQ 
FMESSl  =  0. D  0 
RE  TURN 
END 

DOUBLE  PRECISION  FUNCTION  F S IN P (PN FA , C A, CB , CC ,TH T) 

IMPLICIT  RE AL*8 ( A-h,0-Z  ) 

C  THIS  SUBPROGRAM  COMPUTES  THE  NONLINEAR  SPEED  VOLTAGE  ARISING 
C  FORM  CONSIDERATIONS  OF  SATURATION. 

XI  =  DSIN(THT) 

X2  =  DSIN(THT-2.09  43°5102393dD) 

X3  =  DSIN(THT*2  .0943951 02393D0  ) 

FSIMP=PNFAMCA*X1  4  C3*X2  *CC*X3) 

RE  TURN 
END 

DOUBLE  PRECISION  FUNCTION  F  M  ES  S  2  (  I  L  11 ,  IL  21 ,1  L  33,  I L  4  4  ,  I  L  55, 

1  XLAA, XLAF ,X LAD, 

2  THT,PL1,PL2,PL3) 

IMPLICIT  REAL*MA-M,I,0-Z) 

C  THIS  SUBPROGRAM  COMPUTES  COEFFICIENTS  OF  VOLTAGES  SOURCES  PUT  IN  TO 
C  TAKE  INTO  ACCOUNT  THE  EFFECTS  OF  SATURATION. 

DL  AA=  X  LA  A  *  P  LI *DCOS(THT) 

DL BB=  XLAA*PL1  *DCOS(THT  -  2 . 09 4 395  1 02 3 93 DO ) 

DL  CC  =  XLAA*PLl  *DCO  $(  THT  4  2 . 09 4 39 5  1  C2395 D 0 ) 

DL  A  F  =  X  L  A  F  *  P  L2 
DL  AD=  X  L  AD*  PL3 

FMESS2=  IL11*DLAA  +  IL  2  2*DLBB  4  IL33*DLCC  4  IL44*DLAF  4  IL55*DLAD 
FMESS2  =  O.DO 
RE  TURN 
END 


t 
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APPENDIX  E 


SCR  MODELS  FOR  SPICE  2,  AND  SCEPTRE 

E.l  THE  BIPOLAR  JUNCTION  TRANSISTOR  MODEL 

Development  of  the  SCR  models  begins  with  the  bipolar  junction  transis¬ 
tor  model.  This  BJT  model  is  a  variation  of  the  Ebers-Moll  one  dimensional 
model  in  which  junction  capacitances  are  added  to  account  for  charge  storage 
in  the  depletion  region  and  charge  storage  in  the  base  regions  due  to  diffu¬ 
sion  time  across  the  base  region. 


TRaRICS 


BC 

e 


'BC 


+  CjCo 


1 


11  "  VBC^C 


(81) 


T  BE 

r  trVcs  „  0  .  n 

CBE - e —  e  +  CJEo 


f  ■  W*E 


(82) 


The  equations  describing  the  components  in  the  above  model  are: 


!C  "  *R 


vBE/e  vRF/e 

Ves^  -l)-Irc(eBt  -1) 


lCS' 


*£  =  aR^R  ” 


(84) 


VRf/0  Vnp/ 0 

=  VcS*6  -1)  -  We  BE  -1) 


ES' 


E.2  LIST  OF  PARAMETERS 


ap  =  proportion  of  current  injected  from  emitter  into  the  base  that 
diffuses  to  the  base-collector  junction. 


collector 


FIGURE  23.  MODIFIED  EBERS-MOLL  BIPOLAR  JUNCTION 
TRANSISTOR  MODEL  (INJECTION  MODEL). 
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Tp  =  transit  time  for  carriers  injected  by  the  emitter  into  the 
base-collector  depletion  layer. 

=  base-emitter  saturation  current  for  reverse  biased  base- 
emitter  junction. 

CjEo  =  zero  bias  base-emitter  depletion  layer  capacitance. 

e  =  physical  constant  (at  a  fixed  temperature)  =  0.025  v”^ . 

<J>£  =  emitter  junction  potential  (0.6  -  0.8  v)  -  a  physical  property. 

ctpj  =  proportion  of  current  injected  from  the  collector  into  the  base 
that  diffuses  to  the  base-emitter  junction. 

tp  =  transit  time  for  carriers  injected  by  the  collector  into  the 
base  to  reach  the  base-emitter  depletion  layer  boundary. 

Ip<j  =  base  collector  saturation  current  for  reverse  biased  base- 
collector  junction. 


Cj£0  =  zero  bias  base-collector  depletion  layer  capacitance. 

4>c  =  collector  junction  potential  (0.6  -  0.8v)  -  a  physical  property. 

m  =  the  junction  capacitance  gradient  factor.  For  heavily  doped 

semiconductors  m  =  0.5.  For  lightly  doped  semiconductors  m  =  0.33. 


The  modified  Ebers-Moll  model  as  given  is  the  basic  model  as  applied  to 
SCEPTRE,  it  is  referred  to  as  the  injection  model,  because  the  reference 
currents  Ip  and  Ip  are  the  currents  injected  into  the  base  region  from  the 
collector  and  emitter  respectively. 

SPICE  2  incorporates  a  mathematically  identical  model,  but,  the  SPICE  2 
MODEL  is  referenced  to  the  proportion  of  the  injected  currents  (Ip  ,  Ip)  which 
actually  arrive  at  their  respective  collecting  junctions  (aRIp  =  >  called  Ipp, 
otplp  =  >  called  1^(0 •  For  this  reason,  it  is  referred  to  as  the  charge  trans¬ 
port  model  (meaning  that  the  reference  currents  are  those  actually  transported 
across  the  base  region). 

Utilizing  the  reciprocity  relationship: 


UNCLASSIFIED 

2*  2 


CLEMSON  UNI V  SC  DEPT  OF  ELECTRICAL  AND  COMPUTER  EN--ETC  F/G  10/2 
DYNAMIC  SIMULATION  OF  AIRBORNE  HIGH  POWER  SYSTEMS. <U> 
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AFWAL-TR-80-2115 


E.3  CIRCUIT  SIMPLIFICATION 


The  BJT  model  equations  are  then  rewritten  utilizing  this  relationship 


VBE/e 


VBC/0 


Ic  =  I$(e  -1)  -  Ic/aB(e  ^  -1) 


S  R' 


(85) 


*CC  "  W®!* 


VBC/0 


V6 


IE  =  I$(e  -1)  -  Is/aF(e  DC  -1) 


(86) 


=  !EC  ‘  !CC/aF 


'BC 


t  _B C 
Vs  a  0  4.  r 
—  e  +  CjCo 


1  ‘  VBC^ ^C 


(87) 


"BE 


r  BE 

tfts  e  c 
e  e  jEo 


1  -  VBE/<frE 


m 


(88) 


This  modification  is  advantageous  in  that  one  less  parameter  determina¬ 
tion  is  required  and  that  numerical  methods  are  simplified  due  to  an  in¬ 
crease  in  the  range  of  linearity  of  logarithmic  calculations  based  on  the 
reference  current  I<j  (Reference  33). 


Note.  There  is  no  physical  or  mathematical  modification,  only  a  notation 
modification. 

Still  further  modification  in  form  is  made  for  the  SPICE  2  model  while 
keeping  the  model  mathematically  identical.  The  two  reference  currents 
(current  collected  at  the  emitter  and  current  collected  at  the  collector 
I^q)  are  combined  to  form  a  single  reference  current  (current  collected  total 
Iqj) .  This  I^-j.  is  defined  by 


^BF^0  ®Rf^0 

!CT  =  !CC  ■  !ec  =  h<*  "  e  ) 
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The  combining  of  the  transport  current  into  a  single  term  requires  that  the 
diodes  have  only  base  connections  at  their  anode  points  (NPN  model)  and 
therefore  an  expression  for  base  current  must  be  formed  so  that  the  diode 
currents  may  be  expressed  correctly.  Referring  to  Figure  24. 


(89) 


IB  =  O/otp  -  1 ) Iqc  +  ^aR  ”  ^EC 


(90) 


!b  =  !cc/3f  +  Wbr 


(91) 


The  model  now  appears  as  Figure  25,  the  model  terminal  equations  are  now 


!c  =  !ct  "  !ec/3r 

VBE^9  VBC^6 


(92) 


VBC/8 


=  Is(e  -e  ■"*  )  -  Is/Br  (e  -1) 


*E  =  -ICT  "  !C(/3F 

VBC^8  VBE^0 


(93) 


VBE/0 


=  Is(e  -  e  )  -  Is/eF  le  -1  ) 


XB  =  !CC/BF  +  W^R 

We 


(94) 


V0 


IS/BF  (e  -  1)  +  I$/Br  (e  -  1) 


The  terminal  characteristic  model  used  in  SPICE  2  is  then  given  in  terms  of 

IB  and  Iq  .  This  model  is  shown  in  Figure  26,  with  contact  resistances 

rh>  r.  and  r.  added. 
d  c  e 


E.4  THE  SCR  MODEL  IN  SPICE  2 

The  SCR  model  is  developed  in  SPICE  2  by  connecting  TWO  BJT  models 
(a  PNP  and  a  NPN)  as  shown  in  Figure  27.  In  addition  to  the  two  transistors, 
a  resistor  Rs^un^  is  added  to  the  model.  Since  the  Spice  BJT  model  has  no 
capability  to  simulate  the  reverse  breakdown  of  the  p-n  junctions,  then  a 
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FIGURE  24.  MODIFIED  EBERS-MOLL  BJT  MODEL  (TRANSPORT  MODEL). 


FIGURE  26.  SPICE  2  MODIFIED  EBERS-MOLL  BJT  MODEL. 


diode  DFOR  must  be  added  to  simulate  the  forward  breakover  voltage  of  the 
SCR.  In  addition  DFOR  provides  a  means  of  controlling  the  amount  of  forward 
bias  of  the  center  p-n  junction  of  the  SCR  as  well  as  a  path  of  discharge  of 
the  junction  capacitance  during  that  portion  of  the  SCR  turn-off  transient  in 
which  the  center  junction  is  relieved  if  excess  carriers  in  the  depletion 
region  by  recombination. 

When  using  the  SPICE  2  BJT  AND  DIODE  MODELS,  there  are  a  large  number 
of  parameters  that  may  be  specified  (28  for  the  BJT,  14  for  the  Diodes).  If 
these  parameters  are  not  supplied  via  the  input  deck,  then  SPICE  2  has  built 
in  default  values  to  use  in  order  to  satisfy  the  model  equations.  The  SPICE 2 
SCR  MODEL  equivalent  circuit  is  shown  in  Figure  28  with  applicable  default 
values  applied  to  circuit  elements. 

The  10  circuit  components  of  Figure  28  are  described  by  the  following 
equations  in  which  any  SPICE  2  default  values  used  in  the  model  have  been 
applied  to  the  equations. 


r  =  r 
el  e 


I 


1 


40  V  „  40Vrn  Ic  40V, 


EB,  CB,  AS, 

c  -  h  <e  “  e  )  -  r1*6 

L1  bi  \ 


CBi 


(95) 


-  1)  (96) 


B1  '  Bf 


IS,  40VeBi  M  40VcBl 

Me  1  -  1)  +  g-Me  1 
pr 

1  K1 


-  1) 


(97) 


'BE,”  ~W 


Tf1Is1  40Veb1 


(98) 


'BC,”  4TT 


tr1is1  40Vcb1 


(99) 
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,  •  *«•"; 


40V 


ls  (e 

^diode 


CB, 


1) 


DFOR 


1  +  ^VCB1/VBO^ 


40VBC? 

Ic  (e  2  -  1) 
^diode _ 

1  *  cv/W5 


40Vbe. 

■c2  ■  V6  +1) 


BC. 


■s.  40VBE.  40VU- 

\  *  <e  -  +  \  l*  - 


CBC2  =  CjC2tl_VBC2) 


-0.5 


^shunt '  ^shunt 


From  these  10  equations,  one  develops  the  following  list  of  parameters 
must  be  provided  in  order  for  the  model  to  function. 

1.  r_ 


2. 

3. 

4. 

5. 


Tf, 


6.  Tr 


7.  I, 


diode 


(100) 

(101) 

1102) 

(103) 

(104) 

that 
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8. 


9. 


10. 


The  method  of  determining  these  parameter  values  is  given  in  Reference  24. 


E.5  THE  SCR  MODEL  IN  SCEPTRE 


The  SCR  model  is  developed  in  SCEPTRE  by  connecting  a  PNP  transistor 
and  NPN  transistor  in  the  same  regenerative  feedback  arrangement  as  was  used 
in  the  SPICE  2  model  (Figure  27).  In  the  case  of  SCEPTRE,  however,  the  BJT 
model  has  the  form  shown  in  Figure  23.  Applying  similar  procedures  as  used 
in  developing  the  SPICE  2  SCR  model  leads  to  the  equivalent  circuit  shown  in 
Figure  29.  One  advantage  of  SCEPTRE  is  that  the  ability  to  model  the  de¬ 
pendent  exponential  current  sources  in  the  input  deck  (thereby  not  relying 
on  a  fixed  internal  model  for  the  BJT)  allows  elimination  of  the  diode  used 
in  SPICE  2  to  model  the  SCR  forward  break-over  voltage.  This  SCR  charac¬ 
teristic  is  modeled  by  including  reverse  breakdown  in  the  mathematical  ex¬ 
pression  for  the  base-collector  junction  of  the  transistors. 


It  is  notable  that  a  minor  rearrangement  of  the  form  of  this  model 
[combine  diodes  1^  and  I^  »  combine  capacitors  Cg^  and  Cg^  ,  and  include 
the  junction  capacitor  (defaulted  to  zero)  of  the  gate  to  cathode  junction] 
results  in  the  intrinsic  3-junction  model  of  the  SCR. 


The  SCEPTRE  SCR  model  of  Figure  29  has  13  elements,  but,  since  the 
current  sources  are  defined  as  a  constant,  a,  times  the  corresponding  diode 
current,  then  only  the  equations  for  diode  currents  are  given. 
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SCEPTRE  TWO  BJT  SCR  MODEL. 
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!cs/e  -  i) 

R]  1  +  ^VCB1/VB0-'6 


"BC, 


TR1aR1ICS1  40VCB1 
?0  e 


CBC2  =  CjC2^"VBC2J 
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The  parameters  required  to  implement  the  above 
SPICE  2  parameters  as  follows: 


(107) 

(108) 

(109) 

(110) 

(HI) 

(112) 

(113) 

model  are  obtained  from  the 
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TABLE  3 


CONVERSION  OF  SPICE  2  TO  SCEPTRE 


1. 

S 

r*, 

2. 

'es, 

s 
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* 
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= 
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“R, 

r 
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R1 

6. 
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= 
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7. 

tr 
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* 

tR 

R1 
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9. 
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^ES 
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